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Complex networks, comprised of individual elements that interact with each other through reac- 
tion channels, are ubiquitous across many scientific and engineering disciplines. Examples include 
biochemical, pharmacokinetic, epidemiological, ecological, social, neural, and multi-agent net- 
works. A common approach to modeling such networks is by a master equation that governs the 
dynamic evolution of the joint probability mass function of the underling population process and 
naturally leads to Markovian dynamics for such process. Due however to the nonlinear nature of 
most reactions, the computation and analysis of the resulting stochastic population dynamics is 
a difficult task. This review article provides a coherent and comprehensive coverage of recently 
developed approaches and methods to tackle this problem. After reviewing a general framework 
for modeling Markovian reaction networks and giving specific examples, the authors present nu- 
merical and computational techniques capable of evaluating or approximating the solution of the 
master equation, discuss a recently developed approach for studying the stationary behavior of 
Markovian reaction networks using a potential energy landscape perspective, and provide an intro- 
duction to the emerging theory of thermodynamic analysis of such networks. Three representative 
problems of opinion formation, transcription regulation, and neural network dynamics are used 
as illustrative examples. 
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I. INTRODUCTION 

Complex interaction networks are at the core of many 
problems of scientific and engineering interest, and this 
realization has caused the interdisciplinary stiidy of net- 
works to burgeon over the past decade. Example ap- 
plications include (but are not limited to): chemical 
reaction networks (Heinrich and Schuster, 1996; New- 
man, 2003, 2010), cellular (signaling, transcriptional and 
metabolic) networks (Barabasi and Oltvai, 2004; New- 
man, 2003, 2010), pharmacokinetic networks used to 
study the absorbtion, distribution, metabolism, and elim- 
ination of chemicals and drugs by the human body (Bois 
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et al., 1990), epidemiological (disease-spreading) net- 
works (Hetlicote, 2000; Newman, 2010), ecological net- 
works (Bascompte, 2009, 2010; Black and McKane, 2012; 
Newman, 2003, 2010; Powell and Boland, 2009; Thebault 
and Fontaine, 2010), social networks (Borgatti et al, 
2009; Freeman, 2004; Hill et al., 2010; Masuda et al., 
2010; Newman, 2003, 2010; Weidlich, 2006), neural net- 
works (Benayoun et al., 2010; Newman, 2003, 2010), 
multi-agent networks comprised of intelligent agents that 
observe and act upon each other to achieve a certain ob- 
jective (Xi et al., 2006), and evolutionary game theory 
networks (Szabo and Fath, 2007). 

A common approach to modeling the dynamic behav- 
ior of complex interaction networks is by a master equa- 
tion that governs the time evolution of the joint prob- 
ability mass function of the underling population pro- 
cesses and naturally leads to Markovian dynamics. Due 
however to the nonlinear nature of most interaction net- 
works, computing the exact solution of the master equa- 
tion is not possible in general. As a consequence, the 
analysis of nonlinear Markovian interaction networks is 
a formidable task. Deterministic approximations of the 
master equation have been developed to address this 
problem, but these approximations may fail to predict 
important system behavior (Buice et al., 2010; Gomez- 
Uribe and Verghcsc, 2007; Goutsias, 2007; Leonard and 
Reichl, 1990; McQuarrie et al., 1964; Rao and Arkin, 
2003; Thakur et al., 1978; Vellela and Qian, 2007; Zheng 
and Ross, 1991). For example, deterministic approxi- 
mations cannot predict the emergence of noise-induced 
behavior, a fundamental property of nonlinear interac- 
tion networks with stochastic dynamics (Artyomov et al, 
2007, 2009; Bishop and Qian, 2010; Qian, 2010, 2011; 
Qian et al, 2009; Zhang et al., 2010c). 

The earliest Markovian interaction network model pro- 
posed in the literature seems to be that of Delbriick 
(1940) who developed it to study statistical fluctuations 
in an autocatalytic reaction mechanism of chemical kinet- 
ics. This approach was subsequently adopted by several 
investigators who focused on models of simple reaction 
mechanisms in small systems that exhibit large fluctua- 
tions and developed methods for their analysis (Bartholo- 
may, 1958, 1959, 1962; Darvcy and Staff, 1966; Haken, 
1975; Ishida, 1958, 1964; Kurtz, 1972; McQuarrie, 1963, 
1967; McQuarrie et al., 1964; Nicolis and Prigogine, 1977; 
Schnakenbcrg, 1976; Singer, 1953). Parallel to these de- 
velopments, the pioneering work of N. G. van Kampen 
and D. T. Gillespie provided fundamental analytical and 
computational methods for dealing with stochasticity in 
nonlinear chemical reaction networks through approxi- 
mations of the master equation or Monte Carlo sampling 
(Gillespie, 1976, 1977, 1992, 1996, 2000, 2001; van Kam- 
pen, 1961, 1976, 2007). These methods however were 
largely overlooked by the chemical modeling community 
which, for many decades, concentrated its main effort on 
developing system-based and control-theoretic methods 
for the analysis of chemical reaction networks using de- 
terministic rate equations (Heinrich and Schuster, 1996). 



It turns out that the deterministic approach is theoret- 
ically and computationally miich easier to handle than 
the stochastic approach. Successful application to nu- 
merous chemical modeling and analysis problems is one 
of the main reasons why deterministic approaches have 
garnered wide-spread popularity. 

Strong experimental evidence has recently revealed 
that stochasticity plays a hmdamcntal role in cell regula- 
tion (Blake et al, 2003; Elowitz et al., 2002; Hasty et al., 
2000; Kepler and Elston, 2001; McAdams and Arkin, 
1997; Munsky et al., 2012; Ross et al, 1994; Thattai and 
van Oudenaarden, 2001). This evidence has catalyzed 
a new effort on modeling biochemical reaction networks 
using stochastic (mainly Markovian) approaches, result- 
ing in the development of novel mathematical, compu- 
tational and experimental tools for quantitatively under- 
standing the dynamic interplay between stochastic fluc- 
tuations and system function. In addition to refining 
previously suggested algorithms and developing new nu- 
merical and computational techniques for estimating or 
approximating the solution of the master equation, two 
important and related methodologies are emerging as 
fundamental to the analysis of nonlinear biochemical re- 
action networks. The first is based on a potential energy 
landscape perspective (Ao, 2004; Ao et al., 2007; Han 
and Wang, 2007; Kim and Wang, 2007; Lapidus et al., 
2008; Wang et al, 2010a, 2008, 2010b,c, 2011; Zhou and 
Qian, 2011) and leads to a powerful approach for concep- 
tualizing and quantifying emergent complex behavior in 
nonlinear biochemical reaction networks with stochastic 
dynamics. The second methodology is based on non- 
equilibrium stochastic thermodynamics (Andrieux and 
Gaspard, 2004, 2007; van den Broeck and Esposito, 2010; 
Demirel, 2010; Esposito and van den Broeck, 2010b; Ge, 
2009; Ge and Qian, 2010; Ge et al., 2012; Han and 
Wang, 2008; Jiang et al., 2004; Luo et al, 2002; Mou 
et al, 1986; Puglisi et al, 2010; Qian, 2006, 2009, 2010; 
Rao et al., 2011; Ross, 2008; Ross and Villaverde, 2010; 
Santillan and Qian, 2011; Schlogl, 1980; Schmiedl and 
Seifert, 2007; Schnakenbcrg, 1976; Seifert, 2008; Vellela 
and Qian, 2009; Zhang et al., 2012) and can be effectively 
used to study the macroscopic behavior of Markovian bio- 
chemical reaction networks and, in particular, properties 
related to their self-organization, functional stability, ro- 
bustness and evolutionary behavior (Haken, 1975; Han 
and Wang, 2008; Prigogine, 1978). 

In parallel to the previous developments, substantial 
effort has been independently focused on modeling and 
analyzing stochastic behavior in problems of epidemi- 
ology (Bailey, 1950, 1957, 1963; Bartlett, 1949, 1957, 
1960; Black and McKane, 2010; Chen and Bokka, 2005; 
Haskey, 1954; Hill and Severo, 1969; Jenkinson and Gout- 
sias, 2012; van Kampen, 1973, 1976; Keeling and Ross, 
2008, 2009; Youssef and Scoglio, 2011), ecology (Bartlett, 
1960; Black and McKane, 2012; Datta et al., 2010; Dilao 
and Domingos, 2000; Li et al., 2011), sociology (Haken, 
1975; Weidlich, 1972, 1991, 2006; Wcidhch and Haag, 
1983), and theoretical neuroscience (Benayoun et al.. 
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2010; Bressloff, 2009, 2010; Buicc and Cowan, 2007; Buice 
et ai, 2010; Cowan, 1991; El Boustani and Dcstoxhe, 
2009; Hakcn, 1975; Ohira and Cowan, 1993; Soula and 
Chow, 2007). The main premise underlying this effort is 
the realization that environmental, demographic, behav- 
ioral, and biological factors fluctuate randomly and that 
the resulting stochasticity can cause dramatic deviation 
from what is predicted by deterministic approaches. 

A common theme of most works cited above is the 
representation of stochasticity by a master equation that 
naturally leads to Markovian dynamics. This provides 
a direct mathematical and computational link with the 
techniques developed in stochastic chemical kinetics. As 
a matter of fact, there is a growing consensus among net- 
work researchers in diverse scientific disciplines that most 
mathematical, numerical, and computational tools devel- 
oped for solving problems in stochastic chemical kinetics 
can also be used to solve problems within seemingly dis- 
parate fields of scientific inquiry. It turns out that Marko- 
vian reaction networks provide a unified mathematical 
framework for studying stochastic dynamics on networks 
in a variety of scientific and engineering applications. 

Our main goal in this article is to provide a com- 
prehensive and coherent coverage of recently developed 
approaches and methods to model complex nonlinear 
Markovian reaction networks and analyze their dynamic 
behavior. To achieve this, we first review in Section II a 
general framework for modeling Markovian reaction net- 
works and subsequently discuss specific examples of this 
framework in Section III. In Section IV, we provide a 
comprehensive review of the main numerical and compu- 
tational techniques available for estimating or approxi- 
mating the solution of the master equation. Moreover, 
in Section V, we focus on multiscale methods for ap- 
proximately computing the solution of stiff master equa- 
tions. In addition, we review in Section VI several math- 
ematical facts pertaining to the mesoscopic (probabilis- 
tic) behavior of the master equation. These facts are 
well-known from the theory of Markov processes, but 
we recast them here in the more specific form dictated 
by the framework of Markovian reaction networks. In 
Section VII, we discuss a recently developed approach 
for studying the stationary behavior of Markovian re- 
action networks using a potential energy landscape per- 
spective, whereas, in Section VIII we present an intro- 
duction to the emerging theory of thermodynamic anal- 
ysis of Markovian reaction networks. Finally, we pro- 
vide in Section IX a general outlook of what we be- 
lieve lies ahead in this very fundamental and exciting 
area of research and summarize our conclusions in Sec- 
tion X. To illustrate key concepts, we employ three rep- 
resentative examples dealing with opinion formation in 
social networks, transcriptional control in cell regula- 
tion, and avalanche formation in neural networks. The 
MATLAB software used to implement these examples 
is available on line and can be freely downloaded from 
www.cis.jhu.edu/~goutsias/CSS%201ab/software.html. 



With such a rich and diverse subject matter, the au- 
thors regret that realistic limitations forbid an exhaustive 
treatise on the history and present state of the field. The 
references provided in this review can serve as a start- 
ing point to more in depth or diverse coverage. We sin- 
cerely apologize to the authors whose works do not re- 
ceive recognition, but hope that the listed citations can 
provide a "path of least resistance" to early-stage inves- 
tigators who may feel lost in the vast sea of publications 
available in the area of complex interaction networks. 

II. REACTION NETWORKS 

A. Chemical systems and reaction networks 

Networks of chemical reactions are used extensively to 
model biochemical activity in cells. It turns out that 
many physical and man-made systems of interest to sci- 
ence and engineering can be viewed as special cases of 
chemical reaction networks when it comes to mathemati- 
cal and computational analysis. For this reason, chemical 
reaction networks can serve as archetypal systems when 
studying dynamics on complex networks. 

A chemical reaction system is comprised of a (usually) 
large number of molecular species and chemical reactions. 
A group of molecular species, known as reactants, inter- 
act through a chemical reaction to create a new set of 
molecular species, known as products. In general, we can 
think of a set of chemical reactions as a system that con- 
sists of N molecular species Xi, X2, ■ ■ ■ , that interact 
through M coupled reactions of the form: 

X] ^nmXn ^ ^ ^'nm^n, m € M, (1) 

where A/" := {1, 2, . . . , N} and 7W := {1, 2, . . . , M}. The 
quantities i/„„i > and > are known as the stoi- 
chiometric coefficients of the reactants and products, re- 
spectively. These coefficients tell us how many molecules 
of the n-th species are consumed or produced by the m-th 
reaction. In particular, the notation used in (1) implies 
that occurrence of the m-th reaction changes the molecu- 
lar count of species Xn by 

■^nrn, ■ — ^nm ^nmi wherC SfiYn 

is known as the net stoichiometric coefficient. 

The inter-connectivity between components in a chem- 
ical reaction system can be graphically represented as a 
network (Klamt et al., 2009; Newman, 2010) and, more 
specifically, by means of a directed, weighted, bipartite 
graph. Since molecular species react with each other to 
produce other molecular species, we can refer to this net- 
work in more general terms as a reaction network. 

To illustrate how we can map a chemical reaction sys- 
tem to a network, let us consider the following reactions 
that correspond to a quadratic autocatalator with posi- 
tive feedback (Goutsias, 2007): 
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S 

D + P 
2P 
P + Q 
P 



P 

D + 2P 
P + Q 
2Q 



(2) 



where the last two reactions indicate the degradation 
of molecules P and Q. This chemical reaction system 
is comprised of = 4 molecular species that interact 
through the M = 6 reactions given by (2). We can (ar- 
bitrarily) label the molecular species as = S, X2 — P, 
X3 = D, X4 = Q, and the reactions as 1, 2, . . . , 6. We can 
now represent the system by the network of interactions 
depicted in Fig. 1. This network consists of two types of 
nodes: those representing the molecular species (white 
circles) and those representing the reactions (black cir- 
cles). The directed edges represent interactions between 
molecular species and reactions and, naturally, connect 
only white nodes with black nodes. Edges emanating 
from white nodes and incident to black nodes correspond 
to the reactants associated with a particular reaction, 
whereas, edges emanating from black nodes and incident 
to white nodes correspond to the products of that reac- 
tion. Edges are labeled by their weights, which corre- 
spond to the stoichiometric coefficients associated with 
the molecular species represented by the white nodes 
and the reactions represented by the corresponding black 
nodes. For simplicity, an edge is not labeled when the 
value of the associated stoichiometric coefficient is one. 

An alternative representation of a reaction network is 
by means of the two N x M stoichiometric matrices V 
and V with elements lynm and v'nmi respectively. These 
matrices play a similar role as the adjacency matrix of a 
simple graph (Newman, 2010). For the reaction network 
depicted in Fig. 1, we have that 



V = 



1 

12 110 

1 

1 1 



and 





12 10 

1 

1 2 



It is not difficult to see that, given the two stoichiomet- 
ric matrices V and V, we can uniquely construct the 
chemical reaction system given by (2) and, therefore, the 
network depicted in Fig. 1. Hence, knowledge of the two 
stoichiometric matrices completely specifies the network 
topology. Note that a quick glance of these matrices may 
allow us to make some interesting observations about the 
chemical reaction system under consideration. For exam- 
ple, the fact that all but one of the elements of the first 
row of matrix V are zero indicates that the molecular 
species Xi is a reactant only in one reaction, whereas, 
the fact that the first row of matrix V is zero indicates 
that this species is not produced by any reaction. More- 
over, the last two zero columns of matrix V' indicate that 
reactions 5 and 6 do not result in any products (i.e., they 
act as sink nodes). 




FIG. 1 A directed, weighted, bipartite graphical representa- 
tion of the chemical reaction system given by (2). The molec- 
ular species are represented by the white nodes, whereas, the 
reactions are represented by the black nodes. Edges ema- 
nating from white nodes and incident to black nodes corre- 
spond to the reactants associated with a particular reaction, 
whereas, edges emanating from black nodes and incident to 
white nodes correspond to the products of that reaction. 

Although the mathematical study of the topological 
structure of a reaction network is an important topic of 
research, we will not consider this problem here. More- 
over, we will not consider situations in which the topol- 
ogy of the network varies with time. The reader is re- 
ferred to Newman (2010) and the references therein for 
such topological considerations. Instead, our objective 
is to discuss mathematical methods and computational 
techniques for the modeling and analysis of the dynamic 
behavior of reaction networks. 

B. Stochastic dynamics on reaction networks 

In many reaction networks of interest, the underlying 
reactions may occur at random times. If Zm{t) denotes 
the number of times that the m-tli reaction occurs within 
the time interval [0, i), then {Zm{t),t > 0} will be a ran- 
dom counting process (Ross, 1996). By convention, we 
set Zjn{0) — (i.e., the reaction never occurs before the 
initial time t — 0). We can employ the M x 1 random vec- 
tor Z{t) with elements Zm{t), m — 1,2, ... , M, to charac- 
terize the state of the system at time t > 0. Zm{t) is usu- 
ally referred to as the degree of advancement (DA) of the 
m-th reaction (van Kampen, 2007). For this reason, we 
refer to the multivariate counting process {Z{t),t > 0} 
as the DA process. 

An alternative way to characterize a reaction network 
is by using the x 1 random state vector 



X{t) :=Xo+SZ{t), 



(3) 



for < > 0, where § := V' — V is the net stoichiometric ma- 
trix of the reaction network and Xq is some known value 
of X{t) at time t = 0. Usually, the n-th element Xn{t) 
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of X(t) represents the population number of the n-th 
species present in the system at time t, although this 
may not be true in certain problems (see the examples 
discussed in Sections III-D and III-E). We will be refer- 
ring to the multivariate stochastic process {X{t),t > 0} 
as the population process. For a given initial population 
vector xq, Eq. (3) allows us to uniquely determine the 
random population vector X{t) from the DAs Z{t), pro- 
vided than Z{t) is finite with probability one. 



1. Markovian dynamics 

A large class of reaction networks can be charac- 
terized by Markovian dynamics, in which case we re- 
fer to them as Markovian reaction networks. Marko- 
vian reaction networks are based on the fundamental 
premise that, for a sufficiently small dt, the probabil- 
ity of one reaction to occur within the time interval 
[t, t + dt) is proportional to dt, with proportionality fac- 
tor that depends only on the species population present 
in the system at time t. Specifically, we have that 
Pr [one reaction m occurs within \t,t + dt) \ X{t) = x\ = 
-Km{x)dt + o{dt), for some function Wmix) of the popula- 
tion, known as the propensity function (Gillespie, 2000), 
where o{dt) is a term that goes to zero faster than dt. 
Under these assumptions, {Zm{t),t > 0} is a Markovian 
counting process with intensity T:m{X{t)). In particular, 
the probability p^iz^t) := Pr[Z(i) = z \ Z{0) = 0] as- 
sociated with this process satisfies the following partial 
differential equation (Goutsias, 2005, 2006; Haseltine and 
Rawlings, 2002): 

dpz{z;t) 

dt 

= ^^amiz- em)Pz{z -em;t)- am{z)pz{z;t)Y (4) 
for t > 0, where 

[ TTmixo + ^z), if ^; > 

Q!m(2) ■■= i , . 



and Bm is the m-th column of the MxM identity matrix. 
This equation is initialized by setting pz{z;0) = A{z), 
where A{z) is the Kronecker delta function. It turns out 
that the population process {X{t),t > 0} is a Markov 
process as well with probability px (x; t) := Pi[X{t) = x \ 
X{0) = Xo] that satisfies the following partial differential 
equation (Gillespie, 1992): 

dpx{x]t) 

at 

= '^{'^m{x- Sm)Px{x -Sm;t)- Trm{x)pxix;t)Y (5) 



for t > 0, initialized by px {x; 0) = A(x — Xq), where Sm is 
the m-th column of the net stoichiometric matrix For 
notational simplicity, we hide the dependency of px {x; t) 
on Xq. Most often, Eqs. (4) and (5) are referred to as 
master equations although they are both special cases 
of the well-known forward Kolmogorov equations in the 
theory of Markov processes (van Kampen, 2007) . 

The previous master equations provide a suggestive in- 
terpretation on how the probabilities pz {z; t) and px (x; t) 
evolve as a function of time. For example, Eq. (5) im- 
plies that the probability px (x; t) of the population pro- 
cess X{t) taking value x increases during the time interval 
[t, t+dt) by an amount dtJ^meM '^m{x— 
due to possible transitions from states x — s„i, m e A^, at 
time t, to state x at time t+dt. However, during the same 
time period the probability px{x;t) also decreases by an 
amount dt^^^j^iTjnix)px{x;t) due to possible transi- 
tions from state x at time t to states x + Sm, m e M., 
at time t + dt. Note finally that, in most practical situ- 
ations, the elements of x are limited to being not larger 
than some finite value. As a consequence, we assume that 
Px{x;t) = 'Kjn{x) = 0, when at least one element of x is 
greater than that value. 

2. Hidden Markov models 

Although the DA process uniquely determines the pop- 
ulation process via Eq. (3), the opposite is not true 
in general. This is due to the fact that the matrix 
S^S may not invertible. Invcrtibility of S^S is only 
possible when the nullity of S is zero, in which case 
Z{t) = (S^S)-iS'^[X(t) -Xo] and the DA process can be 
uniquely determined from the population process. There- 
fore, we can consider the DA process to be more in- 
formative in general than the population process. Note 
that, if the solution Pz{z'. t) of the master equation (4) is 
known, then we can calculate the probability mass func- 
tion px{x;t) without having to solve Eq. (5). Since we 
are dealing with discrete random variables, we have that 

Px{x;t) = ^ Pziz;t), (6) 

for t>0, where B{x) {z: x ~ Xq + §2}. 

In many reaction networks, it is much easier to ob- 
serve the population process than the DA process, which 
is usually very difficult or impossible to measure. Thus, 
we can consider the elements of Z{t) as being the hid- 
den state variables of the system under consideration 
and the elements of X{t) as being the observed state 



^ The solution qx(x;t) of Eq. (5) , initialized with an arbi- 
trary probability mass function q{x), is related to the solution 
Px{x;xo, t) of Eq. (5), initialized with A{x — xo), by qxix; t) = 
'Px(x\XQ,t)q{x({) . Therefore, it suffices to only calculate 
px{x;xt:],t), for every Xq such that q{xo) 7^ 0. For this rea- 
son, we focus our discussion on solving Eq. (5) initialized with 
A(x - Xo). 
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variables. If wc choose to model the population process 
by Eq. (3), then we would be using what is known as 
a hidden Markov model (HMM) for our system (Gout- 
sias, 2006) . This opens the possibility of employing well- 
known techniques for the statistical analysis and stochas- 
tic control of HMMs to mathematically and computation- 
ally study stochastic dynamics on reaction networks. 

3. Topological structure and propensity functions 

At a first glance, Eqs. (4) and (5) may give the im- 
pression that the probability distributions Pz{z;t) and 
Px{x;t) of the DA and population processes associated 
with a reaction network do not depend on a detailed 
knowledge of the topological structure of the network. 
This is due to the fact that the previous master equa- 
tions seem to depend only on the difference § = V' — V 
between the stoichiometric matrices V and V and not 
on the individual matrices. This however is not true. It 
turns out that, for all reaction networks encountered in 
practice, the propensity function iTmix) associated with 
the m-th reaction node in the network docs not depend 
on all elements of the state vector x but only on those 
elements associated with the adjacent reactant nodes, as 
specified by the stoichiometric matrix V. In other words, 
the propensity function does not depend on terms involv- 
ing variables on non-adjacent nodes. As a consequence, 
the topological structure of a reaction network directly 
affects its dynamics through this mathematical property 
of the propensity functions. 

III. EXAMPLES 

We now provide a few examples which clearly demon- 
strate that the previously disciisscd general framework 
for reaction networks, based on (1), is sufficiently gen- 
eral to characterize Markovian dynamics on many other 
important networks. Each example is associated with a 
set of "species" that affect each other's population by in- 
teracting through well-defined "reactions." To determine 
the DA and population dynamics, we only need to spec- 
ify the mathematical form of the underlying propensity 
functions from these, the dynamics follow by solving 
Eq. (4) or Eq. (5) for pz{z;t) a.nd px{x;t), respectively. 

A. Biochemical networks 

When dealing with biochemical reactions, wc usually 
assume that the system is well-stirred and in thermal 
equilibrium at fixed volume. It can be shown in this 
case that the probability of a randomly selected combi- 
nation of reactant molecules at time t to react through 
the TO-th reaction during the infinitesimally small time 
interval [t, t+dt) is proportional to dt, with a proportion- 
ality factor Km known as the specific probability rate con- 
stant of the reaction (Gillespie, 1992). As a consequence, 
Pr [one reaction m occurs within [t,t + dt) \ X{t) = x] = 
KmJm{x)dt+o{dt) , where 'ym{x) is the number of distinct 
subsets of molecules that can form a reaction complex at 
time t, given by 



Xfi. 



with [ai > 02] being the Iverson bracket.^ Note that the 
Iverson bracket guarantees that a reaction will proceed 
only if all rcactants arc present in the system. Moreover, 
we use the convention 0! = 1, so (^g") = 1, indicating that 
the rate of a reaction is only determined by the state of 
the reactants. As a consequence, we obtain the following 
propensity functions: 



for m e A4, 



which are said to follow the mass-action law. 

We should note here that certain reactions cannot be 
adequately characterized by propensity functions that 
follow the mass-action law. For example, let us consider 
a reaction Xi + X2 — > ^3 that can occur only when a 
molecule Xi is bomid by at least one molecule X2 at two 
independent binding sites with the same affinity 9. It 
can be shown [e.g., see Dill and Bromberg (2011)] that 
the fraction of molecules Xi bound by X2 is given by 
^a;i/(l -f 6x1). This leads to the following hyperbolic 
propensity function for the reaction: 



ir{xi,X2) 



k6x\X2 
1 + 6x1 



where k is the associated specific probability rate con- 
stant. Clearly, the mathematical form of the propensity 
function of a given reaction depends on the underlying 
molecular mechanism. 

B. Pharmacokinetic networks 

Physiological pharmacokinetic models are used exten- 
sively to study the absorption, distribution, metabolism, 
and elimination of chemicals and drugs by the body of 
animals and humans. As a consequence, they are of 
crucial importance for drug dosing in clinical pharma- 
cology (Hardman and Limbird, 2001). A large class of 
pharmacokinetic models is based on the notion of com- 
partmcntalization (Machcras and Iliadis, 2006). These 
models assume the existence of a central compartment 
(e.g., heart, lungs, brain, etc.), which serves as a site for 
drug administration to peripheral compartments (e.g., 
fat, muscles, central nervous system, and liver). 

To illustrate the connection between pharmacokinetic 
models and Markovian reaction networks, we consider 
here a model for studying the effect of tetrachloroethy- 
lene, a widely used solvent, on carcinogenesis (Bois et al, 



[ai > 02] = 1, if ai > a2, and otherwise. 
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1990). This model assumes a division of the human body 
into the lungs, which serve as the central compartment, 
and four peripheral compartments, namely fat tissue, 
poorly perfused tissue (muscles and skin), richly perfused 
tissue (central nervous system and viscera, except liver), 
and liver. To model this system, we denote by X„ the 
solvent present in the n-th compartment. Then, we can 
represent the system hy N = 5 species interacting by the 
following M = 10 reactions: 



reaction 1 





^ Xi 


reaction 2 




^ X2 


reaction 3 




^ Xi 


reaction 4 




^ X3 


reaction 5 


X3 


^ Xi 


reaction 6 


Xi 


Xi 


reaction 7 


Xi 


^ Xi 


reaction 8 


Xi 




reaction 9 


X5 




reaction 10 


: X5 


^ . 



The underlying reactions model the injection of solvent 
into lung blood (reaction 1), the exchange of one molecule 
of solvent between the lung blood and fat tissue (reac- 
tions 2 & 3), poorly perfused tissue (reactions 4 & 5), 
richly perfused tissue (reactions 6 & 7), and liver tissue 
(reactions 8 & 9), as well as the metabolic clearance of 
the solvent by the liver (reaction 10). 

If we assume that all compartments are homogeneous, 
that the injection of solvent into the lung blood takes 
place at a constant rate ki, and that the probability of 
a randomly selected solvent molecule to move from com- 
partment n to compartment n' within an infinitesimally 
small time interval [t, t + dt) is proportional to dt with 
proportionality constant , then we can model the pre- 
vious pharmacokinetic system as a Markovian reaction 
network with linear mass-action propensity functions 

7ri(a;) = Ki, n2{x) = k^xi, nsix) = K21X2, 

'Ki{x) = K13X1, 'K5{x) = K31X3, ■K6{x) = K14X1, 
'K7{x) = KiiXi, TTsix) = K15X1, ■Kq{x) = K51X5, 

where the n-th element x„ of vector x denotes the 

population of tetrachloroethylene in the n-th compart- 
ment. Moreover, if we assume that tetrachloroethylene 
metabolism in the liver is saturable according to the 
Michaelis-Mentcn relationship of enzyme kinetics (Bois 
et al., 1990), then the propensity function of the last 
reaction will be given by the following nonlinear (hyper- 
bolic) expression (Sanft et al, 2011): 

K +Xr, 

where V, K are two parameters associated with the un- 
derlying metabolic mechanism. 



C. Epidemiological networks 

Epidemiological networks study the spread of infec- 
tious diseases or agents through a population of indi- 
viduals. Although numerous publications can be found 
on the subject, we refer the reader to Newman (2010) 
for an elementary introduction to epidemiological net- 
works. For a mathematical review of deterministic epi- 
demiological models, see Hethcote (2000), whereas, for a 
stochastic modeling approach to epidemiological model- 
ing, see Chen and Bokka (2005). 

To illustrate the connection between epidemiological 
networks and Markovian reaction networks, we consider 
the simplest and most widely used model, known as the 
SIR epidemic model. In this model, an individual in a 
population can be in one of three states with respect to 
a disease: susceptible (S), infected (I), or resistant (R). 
According to this model, there are two types of interac- 
tions that an individual may undergo: (a) if a susceptible 
individual comes into contact with an infectious individ- 
ual, the susceptible person can be infected, and (b) an 
infected individual may become resistant if his immune 
system fights off the infection and confers resistance, or 
if the individual dies by the infection. These interac- 
tions can be modeled by a reaction network comprised 
of = 3 species (S, I, and R) that interact through the 
following M = 2 reactions: 

Xi+ X2 2X2 
X2 ^ X3 , 

where Xi = S, X2 =1 and X3 = R. In this case, 
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We can now assume that the probability of a randomly 
selected susceptible individual at time t to become in- 
fected by a randomly selected infectious individual dur- 
ing an infinitesimally small time interval [t, t + dt) is pro- 
portional to dt, with proportionality factor k\ that does 
not depend on the particular individuals involved. More- 
over, we can assume that the probability of a randomly 
selected infected individual at time t to recover or die 
from the disease during [t, t + dt) is also proportional to 
dt, with proportionality factor K2 that does not depend 
on the particular infected individual. Then, the previ- 
ous interactions lead to a Markovian reaction network 
with mass-action propensity functions given by (Chen 
and Bokka, 2005) 

7ri(a;i,a;2,a;3) = K1X1X2 and ■K2{xi,X2,X3) = K2X2, 

where Xi,X2, X3 are the populations of susceptible, infec- 
tious, and resistant individuals, respectively. 

We can use the previous 3-species/2-reactions motif, 
given by (7), to construct more complex Markovian reac- 
tion networks that model the spread of an infectious dis- 
ease in a population of individuals grouped into classes 
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(e.g., households, work spaces, cities, etc.); see Ben-Zion 
et al. (2010). Wc may group, for example, individuals 
into two classes, those living in Baltimore and Philadel- 
phia, and give each class its own distinct set of vari- 
ables, namely Xi,X2,X3, for susceptible, infected, and 
resistant individuals in Baltimore, as well as X4, X^, Xq, 
for susceptible, infected, and resistant individuals in 
Philadelphia. Each class will be characterized by the 
previous 3-species/2-reactions motif, resulting in the fol- 
lowing four reactions: 



X1+X2 

X4 + X5 
X5 



2X2 
X3 
2X5 
Xe . 



In this case however there is also a flow (by air, road, or 

rail) of individiials between the two different cities, which 
we can model by using the following six reactions: 



Xi 
X4 
X2 

X, 
Xs 
Xe 



X4 
Xi 
X5 
X2 

Xf, 

^3 



The propensity functions associated with these new reac- 
tions will be proportional to the population of the input 
species, with the proportionality factor being the specific 
probability rate constant of an individual traveling from 
one city to the other. In this fashion, we can build com- 
plex Markovian reaction network models for epidemiolog- 
ical dynamics that are more realistic and more predictive 
than traditional deterministic models. 

Likewise, new reactions may be incorporated into the 
epidemiological network to account for additional tran- 
sitions between states. For instance, if we assume that 
a vaccine is available, then we must include the reaction 
Xi in the formulation. Vital dynamics (i.e., births 

and deaths) may also be included in this fashion. For 
example, if infants born at a fixed rate are always sus- 
ceptible, then the reaction — )• Xi must be included in 
the system. Finally, one may consider social networks on 
which epidemiological networks reside. Specifically, age 
stratification in the population (Hethcote, 2000), or the 
scale-free structure of social/sexual networks (Newman, 
2010), may be handled in a manner similar - albeit not 
identical - to the aforementioned geographic considera- 
tions. 

D. Ecological networks 

Ecological networks aim to study interactions among 
organisms living in a particular area as well as between 
these organisms and nonliving physical components of 
the environment, such as air, soil, water, and sunlight. 
The main objective of this type of network is to model 



how mass and energy are transferred from primary pro- 
ducers (or autotrophs), who generate their own energy 
from the sun's rays, up to the apex predators who gather 
their energy and body mass through the consumption of 
prey lower in the food chain. We illustrate here the fact 
that ecological networks can also be modeled as Marko- 
vian reaction networks using a simple example. 

Consider a food web comprised of grass (Xi), rab- 
bits {X2) and wolves (X^), whose net mass at time t 
is given by Xi{t), X2{t) and X3{t), respectively. These 
states can take non-integer values. In particular, Xi (t) = 
X means that, at time t, the mass of grass equals a;-times 
some reference value, and likewise for rabbits and wolves. 
More advanced models may also choose to keep track of 
the number of individuals (Datta et al, 2010). Here how- 
ever we consider a common situation in which the net 
mass of each species is sufficient to describe the system. 

We can assume that changes in mass distribution are 
caused by discrete steps in body size as predators eat 
prey as well as by the mortality that comes with this 
process. In particular, we can model the predation of 
grass by rabbits and rabbits by wolves with the following 
two reactions (Dilao and Domingos, 2000): 

X1+X2 ^ {l + ai)X2 
X2 + X3 ^ (1 + 02)^3, 

where ai , 02 > are constants representing the conver- 
sion factor of mass. Moreover, when rabbits or wolves die 
for reasons other than predation they fertilize the grass. 
We can model this conversion by (Dilao and Domingos, 
2000) 



X2 
X3 



61X1 



where 61 , 62 > are appropriately chosen recycling con- 
stants. As a consequence, the stoichiometric matrices of 
the resulting reaction network, comprised of the N = 3 
species and the M = 4 reactions above, are given by 
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ai -1 -1 
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Under appropriate assumptions, similar to the ones 
made before, the previous interactions lead to a Marko- 
vian reaction network with mass-action propensity func- 
tions given by (Dilao and Domingos, 2000) 

TTi{x) = Ki[xi,X2 > l]a;ia;2, ■!T2{x) = K2[x2,X3 > l]a;2a;3 
773(2;) = K;3[a;2 > l]a;2, 774(2;) = K4[x3 > l]xs, 

where the Iverson brackets are used to make sure that 
the reactions occur only when the net mass of a reactant 
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species is at least as large as the corresponding reference 
value. Here, ki is the specific probability rate constant 
of rabbits eating grass, K2 is the specific probability rate 
constant of wolves eating rabbits, and K3, K4 are the spe- 
cific probability rate constant of natural deaths of rabbits 
and wolves, respectively. 

More complicated ecological reaction network models 
can include geographic considerations, direct competi- 
tion, mutualism, and more complex food chains (Lassig 
et al., 2001; Poweh and Boland, 2009; Thebault and 
Fontaine, 2010). In addition, epidemiological networks 
can be combined with ecological networks to study the 
effects of a disease on a given ecosystem (Auger et al, 
2009). 

E. Social networks 

Recently, interest has emerged in developing mathe- 
matical models for social networks that can be used to 
better understand human behavior. In particular, much 
effort has been devoted to studying dynamics on social 
networks (Antal et al., 2006; Hill et al., 2010; Masuda 
et al, 2010; Moreno et al., 2004; Weidlich, 2006; Zanette 
and Gil, 2006), a problem that has been investigated by 
the physics community many decades ago (Haken, 1975). 
Several models for dealing with dynamic processes on so- 
cial networks are currently available, with many fitting 
nicely into the Markovian reaction framework discussed 
in this review. As an example, we focus on a model of 
opinion formation in social networks, a process that is of 
political, marketing, and general sociological interest. 

The critical behavior of a society moving from a lib- 
eral to a totalitarian political system can be evaluated 
when individuals are endowed with two separate opin- 
ions: a publicly pronounced and a privately held opinion 
for/against the ideology of the ruling party. The pub- 
lic and private opinions of an individual can be different 
when, for example, public dissent against the ruling ide- 
ology is a punishable offence. Along these lines, let us 
consider a fixed homogeneous group of 2L individuals 
who react in the same manner to a given situation. An 
individual simultaneously holds a public and a private 
opinion that each takes values 1/2 or — 1/2 if it is for or 
against the ruling ideology, respectively. Let us denote 
by Xi the net public opinion, which corresponds to the 
sum of the publicly held opinions of all 2L individuals. 
Likewise, let us denote by X2 the net private opinion. We 
are now dealing with N = 2 species interacting through 
the following M = 4 reactions: 



reaction 1 

reaction 2 
reaction 3 
reaction 4 



Xi + X2 ^ 2X1 + X2 

Xi + X2 — > X2 
Xi+ X2^ Xi+ 2X2 
Xi + X2 — > Xi . 



(8) 



The first two reactions model the infiuence of net pri- 
vate opinion X2 on the net public opinion Xi that results 

in a single individual changing her public opinion in sup- 
port of (reaction 1) or against (reaction 2) the ruling 



ideology. In this case, the net private opinion remains 
unchanged, whereas, the net public opinion is increased 
by one in reaction 1 [due to a value change from —1/2 
(against) to 1/2 (for)] and decreased by one in reaction 2 
[due to a value change from 1/2 (for) to —1/2 (against)]. 
Likewise, the subsequent two reactions model the influ- 
ence of net public opinion Xi on the net private opin- 
ion X2 that results in a single individual changing her 
private opinion in support of (reaction 3) or against (reac- 
tion 4) the ruling ideology. These reactions are governed 
by the following propensity functions (Weidlich, 2006): 



TTi{x) = Ki{L — xi) cxp(aia;i -|- a2X2) 
■K2{x) = K\{L + xi) exp(— aia;i — 02X2) 
7r3(a;) = K2{L - X2) exp(a3a;i) 
'K4{x) = K2{L + X2) eTq>{-a3Xi), 



(9) 



where xi, X2 represent the net values of all publicly and 
privately held opinions, respectively, ki,K2 > are two 
specific probability rate constants associated with the 
four reactions, and ai > 0, 02 > 0, as are three model 
parameters. Note that xi and X2 are integer-valued with 
—L < xi,X2 < L, where —L represents total disapproval 
and L represents total approval of the ruling ideology. 

Parameter ai > controls pressure inflicted on public 
opinion due, for example, to oppression of this opinion 
by the ruling party (the value of this parameter is zero 
in the U.S. where free speech is protected, but strictly 
positive in countries where public dissidence has conse- 
quences). On the other hand, parameter 0,2 > controls 
the influence of privately held beliefs on publicly stated 
opinions, whereas, parameter 03 controls how affirmative 
(for 03 > 0) or dissident (for 03 < 0) the private opinion 
is towards the ruling ideology. When the values of ai and 
03 vary, an abrupt change from a liberal to a totalitarian 
political system can be observed (Weidlich, 2006). This 
critical social behavior predicted by the model is reminis- 
cent to the well-known phenomenon of phase transition 
in statistical mechanics and provides a crucial focus of 
study when dealing with opinion spreading in social net- 
works. 

F. Neural networks 

A discussion on reaction networks cannot be complete 

without mentioning biological neural networks. With 100 
billion or more neurons in the human brain connected 
by 100-500 trillion synapses, there is no other reaction 
network that can compete in size and complexity. 

There is a large body of literature surrounding the 
modeling and analysis of biological neural networks. As 
an example, we consider a Markovian reaction model for 
neural networks recently proposed by Benayoun et al. 
(2010) that is intuitive enough for novices in neurobi- 
ology to comprehend and yet rich enough to be a viable 
candidate for understanding many features of this preem- 
inent reaction network. The model consists of L neurons, 
with each neuron being in either a quiescent or an active 
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state. Let X21-1 and X21 denote a quiescent or active 
neuron I, respectively. We can assign the following two 
reactions to the l-th neuron in the network: 

X2I-I + > ^ V111X211 — )• X21 + > l'l'lX2l' 

t^i t^i (10) 

X21 X21-1, 

where Vij measures the synaptic weight between neu- 
rons i and j, with a positive value indicating an exci- 
tatory synapsis and a negative value indicating an in- 
hibitory synapsis. Note that the first reaction models 
transition of the l-th neuron from the quiescent to the 
active state, which is assumed to be influenced by ap- 
propriately weighted active neurons X21', I' ^ I, in the 
network [see Eq. (11) below] that act as "catalysts." On 
the other hand, the second reaction models transition of 
the neuron from the active to the quiescent state, which 
is assumed to occur constitutively. As a consequence, 
we obtain a reaction network with N = 2L species and 
M = 2L reactions. 

We can describe this system by a 2L x 1 state vector 
X with binary-valued 0/1 elements X21-1, X21 indicating 
the state of the l-th neuron (with being quiescent and 
1 being active). Due to the fact that a neuron must be 
either quiescent or active, the state variables must satisfy 
the mass conservation relationships X21-1 + X21 = 1, for 
I = 1,2, . . . ,L. It has been suggested by Benayoun et al. 
(2010) that the probability of the l-th neuron becoming 
active during an infinitesimally small time interval [t, t -\- 
dt), given that the neuron is quiescent at time t, can be 
taken to be X2i-i[(j)i{x) > 0] ta.nh[(pi{x)]dt -\- o{dt), where 
[a > 0] is the Iverson bracket and is the net synaptic 
input to the Z-th neuron, given by 

(/)((a!) = ^ iyi'iX2i' + hi, (11) 
I'^i 

with hi being an external input to the neuron. The 

term X21-1 ensures that the neuron becomes active within 
[t,t dt) only when it is quiescent at time t. As a con- 
sequence, the propensity of the first reaction in (10) will 
be given by 

-^2i-i{x) = X2i-i[(j)i{x) > O]tanh[0;(a;)], (12) 

and therefore depends on the synaptic inputs from neu- 
rons connected to the l-th neuron and any external input 
to that neuron. On the other hand, if we assume that the 
l-th neuron decays from an active to a quiescent state at 
a constant rate 7(, then the propensity of the second re- 
action will be given by 

-^2i{x) = ^1x21, (13) 

where the term X21 ensures that the neuron becomes in- 
active within \t, t -h dt) only when it is active at time t. 

G. Multi-agent networks 

The study of multi-agent networks focuses on systems 
in which many intelligent agents, such as autonomous ve- 



hicles that observe and act upon their environment, inter- 
act with each other to achieve a certain goal. To illustrate 
the fact that multi-agent systems can also be modeled as 
Markovian processes on reaction networks, we consider 
here a system comprised of L autonomous unmanned 
vehicles (AUVs) that can move over a two-dimensional 
bounded rectangular space in a discrete fashion (Xi et al, 
2006). For simplicity, we assume that, at each step, an 
AUV located at a discrete point in space can move 
towards one of four possible directions, namely east to 
point [i -\- west to point {i — north to point 

j + or south to point (i, j — 1) . We want to develop 
a mathematical approach that can be used to describe 
vehicular motion so that the AUVs reach a spatial con- 
figuration X at steady-state with desired probability p{x) 
which assigns high probability over configurations that 
maximize a given design objective and low or zero prob- 
ability over the remaining configurations. The construc- 
tion of such probability can be thought of as an inverse 
problem that can be solved using a statistical mechan- 
ics approach, as the one proposed by Cohn and Kumar 
(2009). 

In the following, we employ two species X21-1 and X21 
whose populations X21-1 and X21 denote the position of 
the Z-th AUV on the two-dimensional rectangular grid. 
For example, if the l-th vehicle is located at point 
on the grid, then X21-1 = i and X21 = j- We can now 
characterize the motion of all AUVs in the multi-agent 
network under consideration hj N = 2L species interact- 
ing through the following M = reactions: 

X2I-I + X2I -\- ^'^{X2l'-l + X2l>) — >■ 

2^2/- 1 + X2I -\- ^ {X2l'-1 -h X2I') 

X2I-I + X21 -\- ^^(X2;'-l + X2i>) — >■ 

X21 + y^(^2;'-i + X21') 
l'^ I 

X2I-I + X2I -\- '^^{X2l>-1 + X2I') — > 

X21-1 + 2X21 -\- ^(Ar2i'_i -h X21') 

i'^ I 

X2I-I + X2I -f '^^{X2l'-l + X2I1) — >■ 

X21-1 + Y.{X2V-1+X2V). (14) 

The first two reactions model one-step motion of the l-th 
AUV towards east /west, whereas, the other two reactions 
model one-step motion towards north/south. Note that, 
when the first reaction occurs, the horizontal position i of 
the l-th AUV is increased by one (transition from A'2(-i 
to 2X21-1), whereas its vertical position j remains un- 
changed (transition from X21 to itself). Moreover, this 
is done by using the positions X2;'_i, X2V , I' 7^ I. of 
the remaining vehicles [see Eq. (16) below], which act as 
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"catalysts" of the reaction. Similar remarks apply for the 
other three reactions as well. 

Let us now define the potential energy V{x) of the 
reaction system being in configuration x at steady-state 

by 



V{x) := < 



-In 



for X gV 



p{xo) ' (15) 
00, otherwise. 



where 2? is a set that contains all permissible vehicle 
configurations (e.g., x should not allow two vehicles to 
occupy the same grid position or positions occupied by 
obstacles, thus avoiding collisions or assignment of ve- 
hicles to grid positions outside the bounded rectangular 
region). Moreover, Xq G V is an appropriately chosen ref- 
erence configuration of zero potential energy. Given that 
X{t) = X, wc can assume that, during the infinitesimally 
small time interval [t, t + dt), the Z-th AUV can move one 
step towards east if two events take place: (a) during 
[t, t + dt), the AUV initiates motion with probability that 
is proportional to dt, with proportionality factor ki, and 
(b) given that the AUV initiates motion during [t, t + dt), 
it moves with probability cxp {—V{x + S4/_3)}, where 
denotes the m-th column of the net stoichiometric matrix 
of the reaction network given by (14). As a consequence, 
the AUV will be moving east with higher probability if 
the motion produces a larger reduction in potential en- 
ergy. Note that parameter m controls the speed of the 
Z-th vehicle, with higher values of ki resulting in faster 
motion. 

By making similar assumptions for vehicle motion to- 
wards the other three directions, the dynamics on the 
reaction network given by (14) will be Markovian with 
propensity functions 



(16) 



for m = 41 - 3,41 - 2,41 - 1,41, I = 1,2,..., L. Note 
that S41-3, S41-2, S4;_i and S41 equal 621-1, -e2i-i, 62;, 
and —e2i, respectively, where is the m-th column of 
the 2L X 2L identity matrix. It turns out that the re- 
sulting master equation governing the population pro- 
cess X{t) has a unique stationary distribution Px(^) •= 
limt^oo Px{x;t), given by the Gibbs distribution 



-V(x) 



where 



-V{x) 



(17) 



(18) 



is the associated partition function. As a consequence 
of Eqs. (15), (17) and (18), we have that Px(^) = P{^)- 
Therefore, the AUVs will asymptotically position them- 
selves in the two-dimensional space at locations x with 
probability p{x), as desired. 



H. Evolutionary game theory 

Game theory deals with mathematical models of con- 
flict and cooperation among intelligent and rational indi- 
viduals. Evolutionary game theory extends the paradigm 
of classical game theory by removing some stringent as- 
sumptions and by naturally incorporating the dynamic 
aspects of learning and experimentation into the prob- 
lem. 

As an example of how evolutionary game theory can 
fit within the current context, suppose that a population 
of L individuals play with each other a game with pos- 
sible strategies Xi,X2, ■ ■ ■ , Xjy. Let Xn{t) be the num- 
ber of individuals playing strategy X„ at time t. Here, 
we consider a simple situation in which each individual 
competes with all other individuals. However, the frame- 
work presented in this paper is also capable of handling 
more general situations, such as those discussed by Sz- 
abo and Path (2007). Given that X{t) =x, let P^ix) be 
the payoff to an individual playing the n-th strategy at 
time t. Based on the current payoff, this individual may 
decide at a random time to follow a new strategy Xn' in 
an attempt to improve his payoff. This can be modeled 
by the following M = N{N — 1) reactions: 



X, 



+ Xn'+^ Xn" 2Xn' + Xn", Tl' ^ 71. 



n"^ n,n' 



Note that, in this case, the number of individuals Xn" 
that follow strategics other than n and n' affect the tran- 
sition of an individual from strategy n to strategy n' [see 
Eq. (19) below] without changing their own strategies 
and, therefore, act as "catalysts." 

There are many alternative propensity functions that 
can be chosen to dictate when players will change their 
strategy, with each corresponding to different learning 
techniques or update rules (Szabo and Path, 2007). A 
common choice however is given by the imitation rule of 
the Moran process (Moran, 1962): 



it{x) 



L T,n"eAf^n"Pn"{x) ' 



(19) 



where k > is a specific probability rate constant detail- 
ing how often individuals choose to update their strate- 
gies. The second term in Eq. (19) is the fraction of in- 
dividuals playing strategy X„, whereas, the third term 
is the fraction of the net payoff paid to individuals who 
play strategy Xn' • These propensity functions have been 
originally developed to model natural selection and ge- 
netic drift in an asexually reproducing population of N 
genetically distinct individuals, where each genotype rep- 
resents a strategy and the payoffs provide measures of 
reproductive fitness. 

I. Petri nets 

Petri nets have been extensively used to describe 
discrete-event distributed systems, a class of systems 
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that are of particular interest in computer science ap- 
plications (Diaz, 2009). A Petri net is a weighted, di- 
rected, bipartite graph, in which the nodes represent 
places and transitions. Places model passive system com- 
ponents, whereas, transitions correspond to events that 
intcr-convert places. Directed arcs join places to tran- 
sitions (connect places that can be converted during a 
transition) and transitions to places (connect a transi- 
tion with the corresponding products). Weights associ- 
ated with arcs indicate the multiplicity of the arc. Each 
place is associated with tokens, indicating the number of 
existing places. Whether or not a transition takes place 
is described by a rule, which may be deterministic or 
stochastic (Diaz, 2009; Haas, 2002), that depends on the 
number of tokens available in the places connecting to the 
transition by incoming arcs. The occurrence of a tran- 
sition results in removing a token from the input places 
and adding a token to the output places of the transition. 

The flow of tokens on a Petri net can be used to model 
the dynamics on a reaction network. As a matter of fact, 
a number of investigators (including Petri himself) have 
already proposed using Petri nets for modeling biochemi- 
cal reaction systems (Chaouiya, 2007; Goss and Peccoud, 
1998; Heiner et ai, 2008; Reddy et al, 1996). This 
approach however is very similar to traditional meth- 
ods for modeling biochemical reaction systems based on 
first-order differential equations or the chemical master 
equation, which have been extensively studied in the lit- 
erature (Gillespie, 1992; Heinrich and Schuster, 1996). 
In particular, Markovian Petri nets are identical to the 
Markovian reaction networks considered in this review, 
with the places playing the role of species and the transi- 
tions representing reactions. It is however important to 
carefully study the theory of stochastic Petri nets (Haas, 
2002), since many results derived in that theory will likely 
prove very useful for the analysis of the Markovian reac- 
tion networks reviewed in this paper. 

IV. SOLVING THE MASTER EQUATION 

Although the algebraic form of the master equa- 
tions (4) and (5) is simple, solving these equations [i.e., 
calculating the probabiUties Pziz; t) and Px{x; t) at each 
time i > 0] is a very difficult task in general. Many meth- 
ods have been proposed in the literature to address this 
problem, which can be grouped into the six general cate- 
gories depicted in Fig. 2. In the following, we discuss the 
most prominent techniques available to date. Whether 
a technique can be applied to a particular problem de- 
pends on the size and complexity of the reaction network 
at hand. 

A. Exact analytical solution 

Deriving exact analytical solutions for Pz{z]t) and 
Px{x;t) is possible only in simple cases [e.g., see Darvey 
and Staff (1966); Gadgil et al. (2005); Gardiner (2010); 
Heuett and Qian (2006); Jahnke and Huisinga (2007); 
Laurenzi (2000); Leonard and Reichl (1990); McQuar- 
rie (1963); McQuarrie et al. (1964); Zhang et al. (2005)]. 
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FIG. 2 Six methods for solving the master equation. Some 
methods can be used to approximate the joint probability 
distributions of the DA and population processes while other 
methods can only be used to approximate marginal distribu- 
tions. Analytical solutions can be obtained only in special 
cases. Numerical solutions are currently limited to small re- 
action networks. Large networks require use of a moment 
approximation scheme or adoption of linear noise approxima- 
tion method as opposed to Monte Carlo sampling. For large 
reaction networks, computing the macroscopic solution may 
be the only feasible choice. This solution however can only 
be trusted at low fluctuation levels. 

For example, an analytical solution for the master equa- 
tion (5) can be derived in the case of a linear reaction net- 
work (i.e., a network with linear propensity functions). 
It has been shown by Gadgil et al. (2005) that, for closed 
linear reaction networks (i.e., linear reaction networks 
with fixed net population), the solution of the master 
equation (5) is a multinomial distribution, provided that 
the initial joint distribution is also multinomial. More- 
over, for open linear reaction networks (i.e., linear re- 
action networks with varying net population), the solu- 
tion of the master equation (5) is a product Poisson dis- 
tribution, provided that the initial joint distribution is 
also product Poisson [see also Heuett and Qian (2006)]. 
These results are special cases of a more general result 
derived by Jahnke and Huisinga (2007) who have shown 
that the probability distribution px (x; t) of the popula- 
tion process in a linear reaction network with initial state 
Xq can be expressed as the convolution of multinomial 
and product Poisson distributions with time-dependent 
parameters that evolve according to well-defined systems 
of first-order linear differential equations [see also (Zhang 
et al, 2005)]. 



13 



B. Numerical solution 

A substantial research effort has recently been focused 
on approximately solving the master equation (5) using 
numerical techniques. Although the methods developed 
so far show promise for addressing this problem, they are 
mostly limited to relatively small reaction networks. For 
this reason, we only provide a brief discussion here. The 
interested reader can find details in the references. 

The master equation (5) can be expressed as a linear 
system of coupled first-order differential equations, given 
by 

^=PpW, (20) 

for t > 0, where p{t) is a x 1 vector that contains 
the nonzero probabilities Px{x;t), x € X, of the popu- 
lation process X{t) and P is a large K x K sparse ma- 
trix whose structure can be inferred directly from the 
master equation. For example, when the columns of 
the stoichiometric matrix § are all different from each 
other, the only nonzero elements of the i-th column of P 
are the M off-diagonal elements, whose values are given 
by TTmixi), and the diagonal element, whose value is given 

by — X]m=i ^™ (•''«)' where M <C -ft' is the number of 
reactions. If we assume that the cardinality K of the 
state-space X is finite, then we can calculate the proba- 
bilities Px{x;t) by solving Eq. (20), in which case 

p(i)=exp(tP)p(0), (21) 

for t > 0. This simple idea has led to a numerical tech- 
nique, proposed by Munsky and Khammash (2006), for 
approximately solving the master equation known as fi- 
nite state projection (FSP). This method requires an ap- 
propriate truncation of the state-space to determine the 
smallest possible set X and development of a compu- 
tationally feasible algorithm for calculating the matrix 
exponential in Eq. (21). 

Although a number of methods are available for com- 
puting matrix exponentials [e.g., see Moler and van 
Loan (2003)], we briefiy discuss here a popular tech- 
nique known as Krylov subspace approximation (KSA) 
method (Sidje, 1998; Sidje and Stewart, 1999). For a 
sufficiently small time step r > 0, this is the best avail- 
able method for approximating the vector p{t + t) = 
exp(TP)p(i), when P is a large and sparse matrix. This 
is done by using a polynomial series expansion of the 
form: 

p{t + t) = copit) + ciTFp{t) + ■■■ + c^„_i(tP)^«- V(t), 

where the coefficients co,c-i, . . . ,Cko-i are estimated 
by minimizing the least-squares error \\p{t + t) — 
p{t + t)\\2- It turns out that the optimal Kq- 
th order polynomial approximation of p{t + t) is a 
point in the ii'o-dimensional Krylov subspace )C{t) = 
span{p(i), TPp(f), . . . , (rP)^o-^p(t)}. This element can 
be approximated by 

p{t + T) := \\p{t)\\2Y{t)exp{TU{t)}ei, 



where Y(t) is a K x Kq matrix whose columns form an 
orthonormal basis for the Krylov subspace JC{t) and M.(t) 
is a Kq X Kq Hessenberg matrix (upper triangular with 
an extra subdiagonal) , both computed by the well-known 
Arnoldi procedure (Sidje and Stewart, 1999). Finally, ei 
is the first column of the Kq x Kq identity matrix. 

The KSA method reduces the problem of calculat- 
ing the exponential of a large and sparse K x K ma- 
trix P to the problem of calculating the exponential of 
the much smaller and dense Kq x Kq matrix H {Kq <C K, 
with Kq = 30-50 being sufficient for many applications). 
Computation of the reduced size problem can be done by 
standard methods, such as a Chebyshev or Fade approx- 
imation (Moler and van Loan, 2003; Sidje, 1998; Sidje 
and Stewart, 1999). Note that we can recursively esti- 
mate the solution p(t) in Eq. (21) at some time tj by 

p(t,)=exp{(t,-i,_i)P}p(t,_i) 

= ||p(i,_i)||2¥(t,_i)exp{(tj -t,_i)H(t,_i)}ei, 

for j = 1,2, . . ., where p(0) = p(0) and = to < < 
t2 < ■ ■ ■ is an increasing sequence of (not necessarily uni- 
formly spaced) time points. These points are selected 
automatically, in conjunction with an appropriately de- 
signed error estimation procedure, to ensure stability and 
accuracy of the overall algorithm (Sidje, 1998). 

Unfortunately, and for most realistic reaction net- 
works, X contains a very large number of states with non- 
negligible probability, thus making the practical imple- 
mentation of FSP difficult. This is a direct consequence 
of the fact that X contains i?i x i?2 x • • • x i?jv distinct el- 
ements, where Rn is an assumed maximum copy number 
of the n-th species. A number of approaches have been 
proposed in the literature to address this problem (Deufl- 
hard et al, 2008; Hegland et al., 2007, 2008; Jahnke and 
Huisinga, 2008; MacNamara et al, 2008; Munsky and 
Khammash, 2007; Peles et al, 2006; Wolf et o/.,'2010; 
Zhang et al, 2010a). Although some approaches per- 
form well, most are limited to small reaction networks. 
It turns out that the most difficult issue associated with 
these methods is solving the resulting system of differen- 
tial equations, which is usually prohibitively large. 

We should point out here that another numerical ap- 
proach has been recently proposed in the literature that 
also attempts to address the previous problem (Jahnke, 
2010; Jahnke and Udrescu, 2010). The method is based 
on representing the probability mass function of the pop- 
ulation process by an appropriately chosen wavelet de- 
composition scheme whose basis elements and the asso- 
ciated wavelet coefficients are being adaptively updated 
in time by solving a much smaller system of linear equa- 
tions. Although preliminary resiilts indicate that the 
method works well, it is not clear at this point whether 
it can be efficiently used to evaluate population proba- 
bilities in reaction networks containing more than a few 
reactions and species. 

The KSA method is based on several approximations, 
whose cumulative effect may appreciably affect its ac- 
curacy, numerical stability and computational efficiency. 
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These drawbacks can be addressed by solving the mas- 
ter equation (4) associated with the DA process, instead 
of Eq. (5). This leads to a recently developed numerical 
technique for solving the master equation known as im- 
plicit Euler (IE) method (Jenkinson and Goutsias, 2012). 
Similarly to the KSA technique, derivation of the IE 
method starts by expressing the master equation (4) as a 
linear system of coupled first-order differential equations, 
given by 



dm 

dt 



Qq{t), 



for t > 0, where q{t) is a Q x 1 vector that contains the 
nonzero probabilities pz{z;t), z G Z, of the DA process 
Z{t) and Q is a large Q x Q sparse matrix whose struc- 
ture can be inferred directly from the master equation 
(each column of Q contains M +1 nonzero elements that 
sum to zero, where M -C Q is the number of reactions). 
Ordering the elements in Z lexicographically results in a 
matrix Q that is lower triangular. As a consequcincc;, and 
for a given time step r > 0, we can use the implicit Eu- 
ler method for solving differential equations (Press et al, 
2007) to estimate q{t) at discrete time points tj :— jr, 

j = 1,2, Thus, given an estimate q{tj-i) of q{tj-\), 

an estimate q{tj) of q{tj) can be obtained by solving the 
following system of linear equations: 

(I-TQ)g(i,)=g(t,-_i), 

where I is the Q x Q identity matrix. It has been shown 
by Jenkinson and Goutsias (2012) that this is possible for 
any value of t and can be efficiently done by a standard 
forward substitution scheme (Press et al., 2007). More- 
over, the resulting method is always stable, producing 
a valid probability vector at each iteration, whereas, its 
accuracy can be controlled by a single parameter, the 
step-size r. Finally, we can use Eq. (6) to obtain an 
estimate p{t) of the probabilities Px{x\ t) from q{t). 

The IE method is computationally superior to KSA 
when the cardinality of the state-space Z is not appre- 
ciably larger than the cardinality of the state-space X . 
This however is not always possible, since the DAs are 
non-decreasing, whereas, the population numbers can ei- 
ther increase or decrease in a way that their values re- 
main within a fixed and bounded domain. As a conse- 
quence, this method can only be used when the number 
of reaction events are sufficiently constrained or remain 
small during a time interval of interest. The IE method 
has been used by Jenkinson and Goutsias (2012) to nu- 
merically approximate the solution of the SIR epidemic 
model discussed in Section III-C with remarkable success 
compared to the KSA method. In this case, the nullity of 
the stoichiometric matrix § is zero and, therefore, there is 
one-to-one correspondence between Z{t) and X{t), which 
implies that the state-spaces Z and X are isomorphic. 

C. Monte Carlo estimation 

Numerical approaches for solving the master equa- 
tion are not practical when the reaction network con- 



tains many reactions and species. In this case, Monte 
Carlo sampling (Liu, 2001) can be used to evaluate 
the statistical behavior of the network. If, by simula- 
tion, we generate L sample trajectories {z^'^\t),t > 0}, 
I = 1,2,...,L, of the DA process {Z{t),t > 0}, then 
we can estimate the dynamics of its moments, such as of 
the means {iJ,z{,m]t) ■= E[^m(i)],i > 0} and covariances 
{czim,m';t) := cov[Zm{t), Zm'it)],t > 0}, by using the 
following Monte Carlo estimators: 



1 ^ 



1=1 



Cz{m,m';t) 



1=1 



Moreover, we can estimate the probability distribution 
Pz{z;t) by using 



1=1 



for t > 0, where A(z) is the Kroncicker delta function. 
Due to the simple relationship between the DA and pop- 
ulation processes given by Eq. (3), we can use similar 
estimators to approximate the dynamic evolution of the 
corresponding population statistics. 

Unfortunately, to obtain sufficiently accurate Monte 
Carlo estimates, we need a large number of sample tra- 
jectories, which is computationally inefficient, especially 
when estimating high-order moments or probability dis- 
tributions.^ This problem can be addressed by devel- 
oping computationally efficient approaches for sampling 
the master equation (4). In the following, we discuss a 
number of methods available in the literature. 

1. Exact sampling 

The simplest way to draw samples from the master 

equation (4) is by using the Gillespie algorithm (Gillespie, 
1976, 1977, 1992). This method can generate a trajectory 
{z{t),t > 0} of the DA process by following two steps. 
First, given that the system is at state z{t) at time t, the 
time t-\-T oi the next reaction to occur can be determined 
by drawing a sample r from the exponential distribution: 



VmeM ) I meM 



am{z{t))}, (22) 



When estimating probability distributions, the issue of efficiently 
sampling low probability events is crucial and becomes the main 
bottleneck for deriving accurate and computationally efficient 
Monte Carlo estimators. 
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rt{m) 



(23) 



for T > 0. Then, which reaction occurs at time i + t can 
be specified by drawing a sample from the probability 
mass function 

Q.rn{z{t)) 

am'{z{t)) ' 

and by increasing the corresponding value of Z by one. 

Unfortunately, this algorithm is computationally de- 
manding, especially when applied to large and highly re- 
active systems, due to the fact that every single reaction 
event must be simulated. Attempts by Gibson and Bruck 
(2000), Cao et al. (2004), and McCoUum et al. (2006) 
to improve the computational efficiency of the Gillespie 
algorithm have produced sampling methods that signif- 
icantly increase computational speed for large reaction 
networks. However, despite these efforts, the previous 
methods are still inefficient, especially when used in con- 
junction with Monte Carlo estimation. For this reason, 
work has focused on developing sampling techniques that 
appreciably reduce computational complexity by trading- 
off accuracy. We discuss some of these methods next. 

2. Langevin approximation 

We can obtain a useful approximation to the master 
equation (4) by assuming that there exists a time step r 
such that, for every time t, two conditions are satis- 
fied: (a) occurrence of reactions within the time inter- 
val [t, t + t) docs not appreciably affect the propensity 
functions am{z{t)), m £ M, and (b) the expected num- 
ber of occurrences of each reaction during [t, t + t) is 
much larger than one. In this case, we can approximate 
the DA process Z{t) by another process Z{t) that satis- 
fies the following equations (Gillespie, 1976, 1977, 1992, 
2000): 



Zm{ij + l)T) = Zm{jT)+amiZ{jT))T+^am{Z{jT))TG^, 

(24) 

for j = 0, 1, . . . , m, <E M, initialized by Z„iiO) = 0, 
for every rri G M. where {Gm,n^ (z A4} are mutually 
independent standard normal random variables that are 
statistically independent of Z. 

We can use Eq. (24) to approximately sample the mas- 
ter equation in an iterative fashion. Starting with zero 
DA values at time zero, we can approximate the DA pro- 
cess at time r by setting Zmir) = am{0)T+^/<XmJp)Tgm\ 

for every m G A^, where gm\ m G A^, are samples inde- 
pendently drawn from the standard normal distribution. 
Then, we can approximate the DA process at time 2t by 
setting Zm{2T) = z^{t) + am{z{T))T + a/q!„(^(t))ts'^\ 

for every m € M., where gm\ m € M, are new samples 
independently drawn from the standard normal distribu- 
tion, and so on. 

Unfortunately, the previous method may result in 
crude approximations of the DA and population pro- 
cesses (Goutsias, 2006). The main culprit is our diffi- 
culty in determining an appropriate time step r so that 



the two required conditions mentioned above are simul- 
taneously satisfied. For example, we may try to reduce t 
so that the propensity functions do not change appre- 
ciably during any time interval [t,t + r), thus satisfying 
the first condition. However, if the reaction network con- 
tains "slow" reactions (a situation that happens often in 
practice), these reactions will occur infrequently during 
the time interval [t,t + r), which will result in violating 
the second condition. Finally, the method may produce 
reaction occurrences within a time interval [jr, {j + 1)t) 
that may result in negative species populations [see also 
the discussion by Melykiiti (2010), pp. 65-71], which may 
not be appropriate in certain types of networks (e.g., in 
biochemical reaction networks). 

It is worthwhile noticing here that, in the limit as 
r — >■ 0+, Eq. (24) converges to the following Langevin 
equations (Gillespie, 1996, 2000): 

dZ„,it) = am{Z{t))dt +\Jam{Z{t))dWra{t), (25) 

for i > 0, TO G M, where {Wm,m G M} are mu- 
tually independent standard Brownian motions whose 
increments {dW„i{t),m G Ai} at time t arc also in- 
dependent of Z{t), which can be used to approximate 
the master equation (4). Note that Eq. (24) provides 
a numerical method for solving the Langevin equations, 
obtained by discretizing Eq. (25) using the well-known 
Euler-Maruyama method (Higham, 2001). For this rea- 
son, the approximation method based on Eq. (24), is 
usually referred to as the Langevin approximation (LA) 
method. Note that, when using Eq. (24), we must make 
sure that r is small enough so that we obtain a good 
approximation to the time-continuous DA process Z{t) 
governed by the Langevin equations. 

Finally, Melykiiti et al. (2010) has recently shown 
that there are many alternative ways to formulate the 
Langevin equations, which result in the same finite- 
dimensional joint probability distribution for the under- 
lying population variables. It turns out that using one 
particular formulation can considerably accelerate im- 
plementation of Monte Carlo estimation. Despite this 
advantage, and for the reasons discussed above, caution 
should be exercised when replacing the master equation 
with the Langevin equations. 

3. Poisson approximation 

The DA process satisfies the following equation (Kurtz, 
1980): 



Zm{t) 



am{Z{t')) dt' 



for t > 0, TO G A^, where P^, to G Al, arc statistically 
independent Poisson random variables with unit rate. As 
a consequence of the Markovian nature of the process, we 
also have that 



Zm{t + T) = Zm{t) + P„ 



J^ZniZit')) dt' 



(26) 
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for t > 0, m, E A4. This result can be used to construct a 
better technique than the LA method for approximately 
sampling the master equation. In particular, we can em- 
ploy a time step r so that occurrence of reactions within 
the time interval [t, t + t) does not appreciably affect 
the propensity functions am{z{t)), m € M. This is the 
first condition required by the LA method, which is com- 
monly referred to as the leap condition. In this case, 
given that Z{t) ^ z{t), the number of occurrences of the 
m-th reaction within the time interval [t, t + t) will ap- 
proximately follow a Poisson distribution with mean and 
variance am{z{t))T. As a consequence, Eq. (26) becomes 

Zmiij + 1)t) = Zmijr) + Pra [am{Z{3T))T^ , (27) 

for j = 0, 1, . . . , me M., initialized by Zm{0) = 0, for 
every m £ M. The resulting method is usually referred 
to as the Poisson approximation (PA) method. 

By using Eq. (27), we expect to obtain accurate sam- 
ples of the DA process, provided that we choose a time 
step r that sufficiently satisfies the leap condition. Hence, 
an important practical problem here is to determine an 
appropriate value for r so that the leaping condition 
is approximately satisfied. We would like this value to 
be as large as possible so that the resulting method is 
appreciably faster than exact sampling using the Gille- 
spie algorithm. Practical considerations however dictate 
that r must not be very large, otherwise the method 
may produce reaction occurrences within a time interval 
[jr, {j + l)r) that may result in negative species popula- 
tions, which is the same problem as the one encountered 
when using the LA method. 

The problem of determining the largest value of r so 
that the leap condition is satisfied has been addressed 
by Gillespie (2001), Gillespie and Petzold (2003), and 
Cao et al. (2006). The latest proccdiirc is accurate, easy 
to code, and results in faster implementation than the 
previous methods. To avoid negative populations, it has 
been suggested by Tian and Burrage (2004) and Chat- 
terjee et al. (2005a, b) to approximate the Poisson distri- 
bution by a binomial distribution. The main rationale 
behind this choice is that the maximum number of oc- 
currences produced by a binomial distribution is always 
bounded and easily controlled by one of the two param- 
eters used to specify the distribution. This however is 
not true for the Poisson distribution, which can produce 
a very large number of occurrences within a small time 
interval (a Poisson random variable takes values between 
and oo) that can falsely result in negative populations. 
Some improvements of the original r-leaping methods 
can be found in Peng et al. (2007) and Pettigrew and 
Resat (2007). 

It turns out that we can still use a Poisson distribu- 
tion for the occurrence of reactions and always guaran- 
tee nonnegative populations. This has been recognized 
by Cao et al. (2005a), who proposed a sampling method 
that is easier to implement than the binomial r-lcaping 
algorithm and is more accurate in general than the orig- 



inal Poisson T-leaping technique. An improved version 
of this approach, which employs a post-leap check to im- 
prove sampling accuracy, has been proposed by Anderson 
(2008). 

Most r-leaping sampling methods available in the lit- 
erature require specification of the mean occurrence of a 
reaction during a leap step. The value used is usually not 
the true mean value and, as a result, a bias is introduced 
that reduces the accuracy and speed of sampling. This 
problem has been addressed in Xu and Cai (2008) by 
specifying an appropriate value for the mean occurrence 
rate obtained directly from the master equation. 

Finally, we refer the reader to Cai and Xu (2007) , Lip- 
shtat (2007), Hellander (2008), Slepoy et al. (2008), Cai 
and Wen (2009), Mjolsness et al. (2009), Ramaswamy 
et al. (2009), and Wu et al. (2011), for alternative simula- 
tion algorithms designed to accelerate exact Monte Carlo 
sampling of the master equation under certain condi- 
tions, as well as to Lu et al. (2004), Anderson (2007), Cai 
(2007), Ramaswamy and Sbalzarini (2011), and Yi et al. 
(2012) for methods dealing with time- varying propensity 
functions and delays. 

4. Weighted sampling 

We can also use Monte Carlo sampling to estimate the 
probability of an event £, whore £ is the collection of 
all trajectories sampled from the master equation that 
satisfy a specific condition of interest (e.g., that the pop- 
ulation Xn (t) of the n-th species exceeds a given thresh- 
old during a time interval [0,to])- If Ti = {zi{t),t > 0}, 
I = 1, 2, . . . , L, are the trajectories obtained by sampling 
the master cqiiation (4), then we can estimate the prob- 
ability of an event £ by employing the following Monte 
Carlo estimator: 



PT[£] = jY,m&£], 



(28) 



1=1 



where [ • ] is the Iverson bracket. 

To produce a sufficiently accurate probability estimate 
when using Eq. (28), we may need to use a prohibitively 
large number of samples, especially when 5 is a rare event 
(i.e., when Pt[£] ^ 1). Rare events are of particular in- 
terest, since they may produce a catastrophic behavior 
in reaction networks, such as the onset of cancer in bio- 
chemical networks or mass population causalities in epi- 
demiological networks. When £ represents a rare event, 
most trajectories sampled from the master equation will 
not be in £ and, therefore, will not contribute to the sum- 
mation in Eq. (28). In this case, we need to appreciably 
increase the value of L in order to accurately estimate 
the probability Pr[f]. 

We can remedy this situation by employing importance 
sampling (Liu, 2001), a classical method for reducing the 
variance of a Monte Carlo estimator and, thus, L. Im- 
portance sampling is based on generating samples drawn 
from a probability distribution which assigns more prob- 
ability mass to trajectories that satisfy the desired condi- 
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tion and less probability mass to the remaining trajecto- 
ries. This approach has been recently employed by Kuwa- 
hara and Mura (2008) to estimate rare event probabili- 
ties in stochastic chemical kinetics and has led to the de- 
velopment and refinement of an innovative approach for 
sampling the master equation, known as weighted sam- 
pling (Daigle Jr. et al., 2011; Gillespie et ai, 2009; Kuwa- 
hara and Mura, 2008; Roh et al, 2010). 

Weighted sampling is based on defining a new set 
{a^(z),m € M} oi propensity functions, given by 
a'jn{z) — XmOim{z), wherc \m,'m G M, are appropriately 
chosen positive constants so that sampling the master 
equation with propensity functions a' produces trajecto- 
ries T/, I = 1,2, ... ,L', which are in £ with high proba- 
bility.^ In this case, the Monte Carlo estimator for the 
probability of event £ to occur will be given by 

1=1 

where wi{T[), I = 1,2, ... ,Z/', are weights that account 
for the bias introduced by sampling the master equation 
with propensity functions a' instead of a. 

To compute the weights, note that a trajectory T/ 
can be specified as T[ = {ti, mi, . . . , r^j , m^, }, where 
nil , ■ ■ ■ ,^K, arc the Ki reactions that occur within the 
time interval of interest [0, t^] and ti, . . . , tki are the time 
steps leading to these reactions. Then, the probability of 
sampling a trajectory T/ — {z'i{t),t G [0,to]} from the 
master equation with propensity functions a' is given by 



i=l 



by virtue of Eqs. (22) and (23), where t; = X^^^i t^. It 
turns out that the weight of each biased trajectory must 
be equal to the ratio of the probability that the trajectory 
was sampled from the master equation with propensity 
functions a to the probability that it was sampled from 
the master equation with propensity functions a'. As a 
consequence. 



Ki 



A. 



E 



1 
Am 



In this way, the weighted sampling algorithm can be used 
to compute intractable rare event probabilities rising the 
previously discussed (exact and approximate) sampling 
techniques, since accurate estimation of such probabili- 
ties usually requires L' <^ L number of sampled trajec- 
tories. 



* Choosing these values requires a great deal of intuition about the 
behavior of the reaction system or advanced algorithmic tech- 
niques, such as those discussed in Daigle Jr. et al. (2011). 



5. Maximum entropy approximation 

As we mentioned before, estimating the probability 
distributions pz{z:t) and px{x;t) by sampling the mas- 
ter equation can be computationally demanding and in 
most cases intractable. Depending on available data, 
the size of the reaction network at hand, and avail- 
able computational resources, it may only be possible 
to accurately estimate the first few moments E[X^(t)], 
k = 1,2, ... ,K, of the population process Xn{t). In 
this case, by invoking the principle of maxim,um, en- 
tropy (MaxEnt), we may be able to approximately de- 
rive an analytical form for the marginal probability dis- 
tribution (.t„; t) := Ea.i....,x„_i,x„+i,...,x„Px(a;;*)- As 
a matter of fact, using MaxEnt to determine an appro- 
priate distribution compatible with given moment infor- 
mation has produced surprisingly good results in many 
diverse scientific disciplines. 

The principle of maximum entropy states that an ap- 
propriate approximation of the true-but-unknown distri- 
bution of Xn{t) is the probability distribution Px{xn',t) 
that maximizes the Shannon entropy 

S{Px;t) = -^Px{x„;t) lnpx{xn;t), 

subject to known information about Xn{t) [e.g., knowl- 
edge of the support of Px{xn',t) and of the moments 
E[X!^(t)], k = 1,2,...,K] (Kapur, 1990; Mead and Pa- 
panicolaou, 1984). This approach is based on a well- 
known principle of scientific objectivity that leads us 
to choose the probability distribution, out of all dis- 
tributions consistent with the given information, which 
maximizes our uncertainly (Shannon entropy) about the 
true distribution. Given the moments E[X^(t)], k = 
1,2, ... ,K, and the fact that x„ is a nonnegative integer- 
valued variable, wc can show that Px{xn\ t) is a univariate 
Gibhs distribution of the form: 



Px{xn;t) 



1 

W) 



exp 




for a;„ > 0, i > 0, where the partition function C^it) is 
defined by 

m :=^exp|-f^Afe(t)^.^|. 

The values of parameters Afc(t), k = 1,2, ... ,K, must be 
chosen so that 

'^x^„Px{xn;t) = ^^x\n;t), 



for t > 0, k — 1,2, . . . , K , where ^^x\'^'it) is the value 
of the fc-th moment of Xn{t) obtained by Monte Carlo 
sampling of the master equation (5) or estimated from 

available data. When only an estimate ]2x\n;t) of the 
mean of the population process Xn{t) is available, the 
MaxEnt approximation oipx{xn; t) is a geometric distri- 
bution, given by 
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Px{Xn;t) 



1 + Mi')(n;t) 



(29) 



for Xn > 0, t > 0. On the other hand, when only es- 
timates of the first two moments of the population pro- 
cess Xn{t) are available, the MaxEnt approximation of 
Px{xn\t) is a quadratic Gibbs distribution, given by 



Px{xn;t)= j ^expj- Ai(i)u- A2(t)u^| 

:[-Xi{t)xn-X2{t)xlY 



I ti>0 

X exp ■ 



for x„ > 0, t > 0. In this case however we need to spec- 
ify the values of the parameters Xi{t) and X2{t) so that 
Px{xn\t) satisfies the underlying constraints imposed by 
knowing the first two moments. Although it is not pos- 
sible to specify these parameters analytically, a num- 
ber of numerical methods, such as the method proposed 
by Bandyopadhyay et al. (2005), can be used to address 
this problem [see also Mohammad- Djafari (1991)]. 

We can extend MaxEnt to deal with multivariate 
marginal probability distributions, such as Px{Xm Xn';t). 
Determining the MaxEnt distribution however becomes 
increasingly difficult as the dimensionality of the proba- 
bility distribution increases (Abramov, 2010). Another 
problem associated with Max;Ent is that the method can 
produce a probability distribution that falsely assigns 
non-negligible probability mass over population values 
that are not stoichiometrically possible [i.e., values that 
do not satisfy Eq. (3)]. We may attempt to address this 
problem by calculating an approximation pz:{z;t) of the 
joint probability distribution pz{z;t) using MaxEnt and 
by then estimating px{x;t) from Eq. (6) by replacing 
Pz{z; t) with pz{z; t). This approach however is only fea- 
sible in the case of small reaction networks that contain 
very few reactions so that estimation of Pz{z; t) by Max- 
Ent is possible. 

6. Stiffness 

In Markovian reaction networks, the firing rates of the 
underlying reactions may vary widely. In this case, most 
computational effort associated with the previous Monte 
Carlo methods will be spent on faithfully simulating the 
firings of fast reactions (i.e., reactions with large propen- 
sity values) , even if simulation of such reactions may not 
be important for determining a particular system behav- 
ior of interest. This leads to stiffness, a serious computa- 
tional problem that results in inefficiently sampling the 
master equation. 

To address stiffness, Rathinam et al. (2003) have pro- 
posed a modified version of the r-leaping method, known 
as im,plicit T-leaping, that allows larger r values to be 
used when applied to stiff reaction networks than the 
original r-leaping method (which must use a small step- 
size in this case). Subsequently, Cao et al. (2007) pro- 
posed an adaptive method that identifies stiffness at each 



simulation step and automatically chooses between the 
standard and implicit T-leaping methods. Moreover, it 
provides an appropriate value for r to be used during 
each iteration. 

Although both approaches can appreciably decrease 
simulation time when compared to the standard r- 
leaping method, they may excessively damp fluctuations. 
As a consequence, these methods may underestimate the 
population variances and produce stochastic dynamics 
that evolve tightly around their means. This may not 
be accurate, especially when dominant stochastic fluc- 
tuations are present in the system. To ameliorate this 
problem, Rathinam et al. (2003) have proposed a strat- 
egy that attempts to restore overly damped fluctuations. 
It is not clear however whether this strategy performs 
well when used in more complex reaction networks than 
the simple networks considered by the investigators. 

Another problem associated with the previous tech- 
niques is their difficulty in effectively dealing with very 
small populations of species involved in very fast reac- 
tions. Moreover, these methods may lead to non-integer 
and possibly negative populations, which may not be 
physically meaningful in cc;rtain types of networks (e.g., 
in biochemical reaction networks). Although a stoichio- 
metrically consistent rounding step has been proposed 
by Rathinam et al. (2003) to remedy the last problem, 
it has been observed that rounding may seriously impair 
the performance of the resulting algorithm (Rathinam 
and El Samad, 2007). To address this issue, Rathinam 
and El Samad (2007) have proposed a method based on 
decomposing a Markovian reaction network into "motifs" 
and on constructing appropriate approximations for each 
individual motif. This method however is cumbersome 
and difficult to use, since its effectiveness relies heavily 
on identifying appropriate "motifs," a task that may not 
be possible in large reaction networks. 

Finally, a "partitioned leaping" approach has been pro- 
posed by Harris and Clancy (2006) [see also Harris et al. 
(2009)] that is closely related to r-leaping. At each step, 
the algorithm uses the expected number of firings of a 
reaction within a calculated step-size r to classify the 
reaction into four distinct categories: very slow, slow, 
medium, and fast. Based on this classification, the sim- 
ulation of very slow reactions proceeds using exact sam- 
pling. On the other hand, slow and medium reactions 
are simulated using the Poisson and Langevin approxi- 
mations, respectively. Finally, the number of firings of a 
fast reaction is specified deterministically by multiplying 
the propensity function of the reaction with r. 

"Partitioned leaping" is an attractive idea for 
speeding-up Monte Carlo sampling of the master equa- 
tion. Its accuracy however depends on correctly clas- 
sifying the reactions, which may not be always possi- 
ble. Although the simple examples provided by Harris 
and Clancy (2006) and Harris et al. (2009) demonstrate 
its effectiveness for reducing computations while preserv- 
ing accuracy, it is not clear how the method will per- 
form when applied on larger and more complex networks 
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with widely disparate reaction rates and how robust the 
method is to possible misclassification of reactions. 

Stiffness in Markovian reaction networks is an impor- 
tant problem, directly related to the practical use of these 
networks. Unfortunately, no sufficient solution has been 
proposed to date a nd more research is needed to satis- 
factorily address this problem. We will revisit stiffness in 
Section V when we discuss multiscale approximations to 
the master equation. 

D. Moment approximation 

When a reaction network contains many species and 
reactions, it may not be possible to accurately estimate, 
in a reasonable time, the statistical behavior of the DA 
and population processes by Monte Carlo sampling. In 
this case, we may try an alternative technique known as 
moment closure. This method can be derived by setting 



ZUt)=^^zim■,t)+Wrn{t), 



(30) 



ioT t > 0, m & M, where iizijn;t) is the mean of Z^{t) 
and Wm{t) := Zm{t) — ^z{m;t). Note that Wm{t) is 
additive zero-mean noise that quantifies fluctuations of 
the DA process around its mean. By using the mas- 
ter equation (4), we can show that the means Hzi'm^t) 
and covariances Cz{m,m';t) := cav[Zm{t), Zm'{t)] = 
^[Wm{t)Wm' (t)] of the DA process satisfy the following 
system of first-order differential equations: 



diiz{rn;t) 
dt 

for t > 0, m e M., and 



(31) 



dcz{m,m';t) 
dt 



E[a„(Z(t))]A(m-m') 

+ E[[Z„(t) - nz{m;t)] am'{Z{t))] 

+ E[[Zm'{t) - Hz{m';t)] a„^(Z(^))], (32) 

for t > 0, m, m' e M., where A is the Kronecker 
delta. Note that the derivatives diiz{m\t) / dt and 
dcz{rn,m' ]t) / dt always exist at finite times, since the 
master equation (4) is valid only when the probability 
mass function of the DA process is a continuous function 
of t. This implies that the means and covariances will 
also be continuous in t and, thus, differentiable. 

In general, we cannot derive an exact solution to the 
previous equations, unless we employ some approxima- 
tion. Since Z{t) = Hzit) + W(t), we can use the follow- 
ing Taylor series expansion of am{Z{t)) around the mean 
value fiz{t)- 



Q„(Z(t)) ~ a™(/iz(i)) + ^™,mi(Mz(i)) W™i(i) 



+ II Vmi,m2(Mz(i)) W^mi(*)W"m2(*) 

"'~5 X/ X/ H ''■'".'"1 



m2,TO3 il^zit)) 

mi£M m2£A1 msGM 

xWmdt)Wm,{t)Wm,{t), (33) 

where /im.mi (A*z(i)) is the first-order derivative of Q.„i{z) 
at the mean value fizit), whereas hm,mi.m2{lJ-z{t)) 
and /iTO,mi,m2,m3(Mz(i)) are the second- and third-order 
derivatives, respectively.^ Then, Eqs. (31)-(33) imply 
that 



diJ,z{m;t) 



dt 



Q!m(Mz(i)) 



^ X] X] '^m,mi,m2{f^z{t)) Cz{mi,m2;t) 



c X! H H ^"i,mi, 



mi&M m-i&M ms&M 



X Cz{mi,m2,ms;t), 



(34) 



for f > 0, m e A4, and 

dcz{m, m'; t) 



dt 



am{fJ'z{t))A{m - m') 



+ X] [f^rn',mi{lJ'z{t)) Cz{m,mi;t) 

+ hm,miilJ'zit)) Cz{m\mi]t) 
"•■^[H H f^m,mi,m2{f^z{t))Cz{mi,m2;t) 

X A(m - m ) 

,mi,m2 

{fJ.z{t)) Cz{m,mi,m2;t) 

mieMm2GM 

~t~ ^m,mi,m2 iHzit)) Cz{m ,mi,m2;t) 
(Mz(t)) 

mi&Mm2&Mm3&M 

X Cz{mi,m2,ms;t) A(m — m') 



"'"5 X/ H X/ \^''m',mi,m2,m3{l^z{t)) 

miGMm2GMm3GM 

X Cz[m,mi,m2,m3;t) 

~i~ ^m,mi,m2,m3 

{fjiz{t))cz{m',mi,m2,ms;t)\, (35) 



For simplicity, wc assume here that the propensity functions are 
sufficiently smooth so that the derivatives of order > 4 are all 
negligible. This condition is satisfied in many cases of interest. 
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for t > 0, TO, to' € A4, where Cz(toi, TO2, TO3; t) and 
Czijni, TO,2, TO3, TO4; are the third- and fourth-order cen- 
tral moments of Z{t), respectively. Eqs. (34) and (35) 
show that the mean and covariance dynamics of the DA 
process Z{t) are in general governed by a system of cou- 
pled first-order differential equations driven by the third- 
and fourth-order central moments. Moreover, the depen- 
dency of these equations on propensity function deriva- 
tives tells us how the mean and fluctuation dynamics arc 
affected by the presence of nonlinearities in the propen- 
sity functions. 

Knowledge of the mean and covariance dynamics 
{fj.z{t), Czit), t > 0} oi the DA process allows us to 
directly calculate the mean and covariance dynamics of 
the population process by fJ.x{t) = Xq + S/iz(i) and 
Cx{t) — SCz(i)S^, respectively. These relationships 
can be used to derive differential equations similar to 
Eqs. (34) and (35) that govern the mean and covariance 
dynamics associated with the population process. 

If the propensity functions arc all linear, then the 
means can be calculated independently from the covari- 
ances and, in general, the equations for the fc-th order 
moments will decouple from all moments of order greater 
than k (Engblom, 2006). In all other cases however, eval- 
uation of the mean and covariance dynamics using the 
previous differential equations requires calculating the 
dynamics of at least third-order central moments. In the- 
ory, these dynamics can be evaluatcid by differential equa- 
tions similar to the ones above, which require evaluation 
of higher-order central moments. However, calculating 
high-order moment dynamics is a formidable task, espe- 
cially when dealing with large reaction networks. Note 
that a reaction network comprised of M reactions re- 
quires M equations for the DA means, M{M + l)/2 
equations for the DA covariances and, in general, 0{M'^) 
equations for the fc-th-order DA moments. As a conse- 
quence, including differential equations governing the dy- 
namics of third- and higher-order central moments of the 
DA process may make sense only when dealing with very 
small reaction networks. Therefore, a practical treatment 
of reaction networks by calculating moments is most of- 
ten limited to evaluating only the mean and covariance 
dynamics. 

When the propensity functions are nonlinear, the mo- 
ment equations always form an infinite hierarchy, with 
lower order moments depending on higher order mo- 
ments, indicating that exact solutions cannot be obtained 
in practice. To address this problem, we can replace the 
moments at some stage of the hierarchy with appropri- 
ately chosen functions of lower-order moments (Gillespie, 
2009; Keeling, 2000; Krishnarajah et ai, 2005; Murrell 
et al, 2004; Whittle, 1957). This approach results in a 
moment closure scheme that produces a self-contained 
system of differential equations whose solution provides 
approximate values for the moments of the DA and pop- 
ulations processes. The resulting method is usually re- 
ferred to as the moment approximation (MA) method. 



A method to construct appropriate functions for mo- 
ment closure is based on making an ansatz for the joint 
probability distributions of the DA or population process. 
For example, we may assume a joint probability distri- 
bution for the DA process that can be uniquely specified 
from the means and covariances. This distribution will 
then impose functional relationships between the third- 
order central moments and the means and covariances of 
Z{t), which can be used to close the system of Eqs. (34) 
and (35) . Evaluation of the mean and covariance dynam- 
ics will now require solving a coupled system of first-order 
differential equations comprised of M equations for the 
means and M{M -I- l)/2 equations for the covariances. 

To illustrate this method, let us consider a reac- 
tion network whose propensity functions are at most 
quadratic. If we make the ansatz that the probability 
distribution of the DA process Z{t) is approximately nor- 
mal, then the third-order central moment will be zero 
and, in this case, Eqs. (34) and (35) will be exact, result- 
ing in 



dfj.z{m\t) 
Jt 



= am(Mz(0) 

for t > 0, TO G M, and 
dcz{m, to'; t) 



Cz{mi,m2;t), (36) 



dt 
1 



[a™(Mz(t)) 



+ 9 E E ^rn,mi,m2Cz{mi,m2;t) A{m-m') 



E [^m',mi(Mz(i)) Cz(TO,TOi;t) 



+ hm,mi{(J'z{t))cz{m',m-i;t) , 



(37) 



for t > 0, TO, to' e A^, where the second-order derivatives 



hr. 



of the propensity functions do not depend on 



the means. More details on this normal MA scheme can 
be found in Goutsias (2007), Gomez-Uribe and Verghese 
(2007), Ferm et al. (2008), El Boustani and Destexhc 
(2009), Lee et al. (2009), Ullah and Wolkenhauer (2009), 
Lafuerza and Toral (2010), and Milner et al. (2011). 

In addition to the normal distribution, a number of 
alternative approximating distributions have been sug- 
gested in the literature, such as log-normal (Keeling, 
2000; Krishnarajah et al., 2005; Nasell, 2003a), Pois- 
son (Nasell, 2003a) and near-Poisson (Buice et al., 2010), 
binomial (Kiss and Simon, 2012; Nasell, 2003a), beta bi- 
nomial (Krishnarajah et al., 2005), and mixtures of distri- 
butions (Krishnarajah et al., 2005; Krishnarajan et al., 
2007). Using the log-normal distribution has some ad- 
vantages over using the normal distribution (Keeling, 
2000; Krishnarajah et al., 2005; Nasell, 2003a). In par- 
ticular, the log-normal distribution has nonnegative sup- 
port and exhibits nonzero skewness, two properties that 
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are important in the context of certain nonlinear reaction 
networks, such as biochemical reaction networks. On the 
other hand, the beta binomial distribution is a discrete 
distribution with a flexible shape that, in some cases, 
can capture the dynamic evolution of the true probability 
distribution of the population process better than other 
distributions (Krishnarajah et al., 2005). Its use however 
is limited to closed reaction networks (i.e., reaction net- 
works with fixed total population) that contain only two 
species. Finally, using mixture distributions for deriving 
a moment closure scheme shows great promise but has 
only been employed in few limited cases (Krishnarajah 
et al, 2005; Krishnarajan et at, 2007). 

We should note here that it may be difficult to spec- 
ify an appropriate probability distribution for the pop- 
ulation process, since this distribution must assign zero 
probability to stoichiometrically impossible populations 
[i.e., populations that do not satisfy Eq. (3)]. On the 
other hand, it may be easier to specify a probability dis- 
tribution for the DA process, since this process is usually 
confined within a well-defined subset of the positive or- 
thant of the multidimensional DA state-space. Note also 
that approximating the moments of the DA process by 
employing continuous distributions, such as log-normal, 
may be difficult to justify due to the discrete nature of 
this process (see however the following siibscction for a 
case in which this may be possible). But more impor- 
tantly, assuming a specific form for the probability dis- 
tributions of the DA process may lead to serious prob- 
lems, since the differential equations derived from the 
master equation will not be consistent with the moment 
structure imposed by the assumed distribution, unless 
the solution to the master equation coincides with that 
distribution. To ameliorate this problem, note that we 
do not necessarily need to specify the exact probabil- 
ity distribution for the DA process in order to close the 
system of moment equations. For example, we can use 
Eqs. (34) and (35) and assume that, for every t > 0, the 
third- and fourth-order central moments are related to 
the first- and second-order central moments by the same 
formulas as the ones associated with a multivariate nor- 
mal or log-normal distribution (Isserlis, 1918; Lafuerza 
and Toral, 2010; Skoulakis, 2008). This is a weaker ansatz 
than assuming that the probability distribution of the 
DA process is multivariate normal or log-normal, which 
may work well in certain circumstances. As a matter of 
fact, it has been shown by Singh and Hespanha (2007, 
2011) that, under certain conditions (at most quadratic 
propensity functions and a moment closure formula that 
has a particular separable form), this assumption is a 
consequence of matching time derivatives of the exact 
(not closed) moment equations at the initial time t = 
with that of the approximate (closed) moment equations. 

Although, in some problems, the previous strategy 
may lead to a sufficiently accurate estimation of the low- 
order moments of the DA and population processes, it 
cannot provide an analytical expression for the proba- 
bility distributions of these processes. We can however 



address this problem by using the MaxEnt approach dis- 
cussed in Section IV-C-5. For example, if the propen- 
sity functions of the reaction network at hand arc at 
most quadratic and if we employ a log-normal-based mo- 
ment closure scheme, then the MaxEnt approximation of 
the probability distribution Pz{zm;t) of the DA process 
(t) associated with the m-th reaction will be given by 
the following Gibbs distribution: 

Pz{zm;t)=(^ J2 exp {- Ai(t)u - X2{t)u^ - A3(i)w^}) ^ 

>[-Xi{t)z^-X2it)zl-X3{t)zi}, 



u>0 

X exp ■ 



where the coefficients Xi{t), X2{t), and Xs{t) must be 
determined so that 



ZmPz{Zm;t) = Hz{m;t) 
X] z'^Pz{zm;t) = Cz{m,m;t) + iil{m;t) 



Zr,.>0 



Cz{m,m;t) + ii%{m;t) 



n 3 



2m>0 

where the last constraint is due to the fact that 

E[Z,{t)Z2{t)Zs{t)] 

^ E[Zi{t)Z2(tmZi{t)Zs{t)MZ2it)Z-i{t)] 

E[z,(tmz2{tmzs{t)] 



(38) 



for the adopted log-normal-based moment closure 
scheme. In this case, the mean and covariance dynam- 
ics can be calculated from Eqs. (34) and (35), by setting 
the fourth-order derivatives h 



mi,m2,m3,m4 



equal to zero 

and by using the relation in Eq. (38) for the third-order 
moments. 

Another approach to deal with the problem of moment 
closure is to replace the true (but unknown) values of 
the higher-order moments required by the MA method 
[such as the third- and fourth-order central moments in 
Eqs. (34) and (35)] with estimated values derived from 
available data or from sampling the master equation us- 
ing Monte Carlo (Chevalier and El-Samad, 2011; Ruess 
et al., 2011). A technique proposed by Ruess et al. 
(2011) employs a small number of Monte Carlo samples 
to obtain crude estimates of the higher-order moments. 
The resulting estimates are interpreted as noisy mea- 
surements of the true moment values and an extended 
Kalman filtering approach is then used to obtain more 
accurate estimates of these values. On the other hand, 
a strategy proposed by Chevalier and El-Samad (2011) 
replaces the unknown moment values with estimated val- 
ues obtained from available data. Both approaches may 
work well when the estimation error is small. However, 
large errors may lead to erroneous calculations due to 
potential amplification of these errors when solving the 
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differential equations that govern the dynamic cvohition 
of the lower-order moments [such as Eqs. (34) and (35)]. 
To address this problem, a large number of Monte Carlo 
samples may be needed when estimating the higher-order 
moments, which can appreciably decrease the computa- 
tional efficiency of the first method. On the other hand, 
the second approach requires a large amount of data to be 
available for reliable estimation of moments and, there- 
fore, it is limited to a small number of problems in which 
this may be possible. 

E. Linear noise approximation 

In certain circumstances, the joint probability distribu- 
tions of the DA and population processes can be well ap- 
proximated by multivariate normal distributions. To see 
why this is true, we will assume the existence of a system 
parameter fl that measures the relative size of stochas- 
tic fiuctuations in a Markovian reaction network, such 
that fluctuations are small for large fi. This is motivated 
by the fact that, in chemical reaction systems, stochastic 
fiuctuations gradually diminish as the system approaches 
the thermodynamic limit at which the population of each 
species and the system volume approach infinity in a way 
that the concentrations remain the same.^ 

It is intuitive to expect that the probability of a re- 
action to occur within the infinitesimally small time in- 
terval [t,t + dt) depends on the "density" x{t)/^ of the 
population process at time t and that this probability 
does not change when Q varies as long as the population 
densities remain fixed (van Kampen, 2007) . This implies 
that the propensity functions Wm must satisfy 77^ (a;; SI) = 
TT„i{x/n), where tt„i docs not depend on flJ To be more 
general, we may also add a term Sl^^TT'^{x/fl), in which 
case we would like iTmix; fl) = 7r„(x/0) + f7^^7r^(x/0).^ 
Moreover, we can assume that ??„( • ) and 7f^( • ) are an- 
alytic. Finally, we may allow an arbitrary positive factor 
/(O), such that 

TTmix; O) = /(fi) [nm{x/n) + Q-^w'^ix/^)] . (39) 

As a consequence, we can assume the following scaling 
law for the propensity functions of the DA process: 

am{z; Q) = f{fl) [am{z/n) + n-^a'„{z/n)] , (40) 

for m e M, where ajn{z/V,) := Tr„(xo/fi -I- Sz/Q) and 
a'^iz/n) ■.= n'^{xo/n + nz/n). 

To proceed, we can make the following ansatz: 

Zm{t;n) = Crn{t) + ^Sm{t), (41) 



® We simply denote the thermodynamic limit as Q — ^ oo. 

^ When necessary, we explicitly denote the dependance of various 

quantities on Q. 

® As a matter of fact, wc can also add terms of 0{fl^^), 0{fl~^), 
etc., if necessary. These terms can be easily accommodated in 
the formulation. However, it is not necessary to do that here. 



for t > 0, m G Ai, where Zm{t;il.) is the "density" 
Zm{t;i^)/Sl, of the DA process, Em{t) is a noise com- 
ponent that quantifies the fluctuations associated with 
the DA process, and Cm(i) is a deterministic process that 
satisfies: 



dCmjt) 
dt 



am{C{t)), 



(42) 



for t > and m G M., initialized with Cm(0) — 0. For 
each SI, Eq. (41) decomposes the random DA density 
Zm{t;Sl) into a deterministic component Cm(t) and an 
additive noise component Sm(t). Clearly, this equation 
is based on the premise that the fiuctuations diminish to 
zero as fast as f2~^/^. In contrast to Eq. (30), which is 
exact, Eq. (41) must be justified. This can be done by a 
central limit theorem for the behavior of the probability 
density function of the DA density process Z{t;Sl), as 
SI — > oo, similar to that shown by Kurtz (1971, 1972) for 
the case of biochemical reaction networks. 

By using Eqs. (40) (42) and the fi-expansion method 
of van Kampen, it can be shown that, for a sufficiently 
large SI, the dynamic evolution of the probability den- 
sity function ps{^;t) of the noise vector H(f) is approx- 
imately governed by the following linear Fokker-Planck 
equation (van Kampen, 1961, 1976, 2007): 



dt 



meAi 



■E E 

meAi m'eAi 



daraim) d[U'P^{i:t)] 



for t > Q, initialized with J3h(^;0) = (5(^), where 5{-) is 
the Dirac delta function. In this case, H(t) will approx- 
imately be a normal random vector with zero mean and 
correlation matrix Ch {t) that satisfies the following Lya- 
punov matrix differential equation: 



d£^{t) 
dt 



A{t) + G{t)C^{t) + C^{t)G^{t), (43) 



for t > 0, initialized with Ch(0) = 0. In this equation, 
A(t) and G{t) arc two M x M matrices with elements 

am,m'{t) = amiCit)) A(m - m') 



'(t) 



dCm' 



respectively, where A is the Kronecker delta function. 
As a consequence, and for a sufficiently large fl, we can 
approximate the probability distribution pz(z'.t) of the 
DA density process by a multivariate normal probability 
density function with mean ^(t) and covariance matrix 
Cs{t)/S}. Due to Eq. (3), this also allows us to approx- 
imate the probability distribution {x; t) of the pop- 
ulation density process X{t;Sl) := X{t)/Sl by a multi- 
variate normal probability density function with mean 
Xq/SI -\- S^(t) and covariance matrix 'B'Cs{t)'B^ ■ As a 
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consequence, and since Z{t) = QZ{t), we can also ap- 
proximate the probability distribution Pz(z; t) of the DA 
process with a multivariate normal distribution, with 
mean il^it) and covariance matrix QCs{t), whereas, we 
can approximate the probability distribution px{x;t) of 
the population process with a multivariate normal dis- 
tribution with mean Xq + i^S^{t) and covariance matrix 

Because the fluctuations in the reaction network arc 
governed by the linear "signal-plus-noise" ansatz given 
by Eq. (41), the previous method is known as linear noise 
approximation (LNA). Its use requires specification of 
an appropriate fluctuation size parameter fl, such that 
Eq. (41) is satisfied, and a sufficiently large value for 
this parameter so that the method produces a reason- 
able approximation of the two probability distributions 
Pz{z;t) and Px{x;t). Implementation of the method re- 
quires one to separately solve the system of M first-order 
differential equations (42) and the system of M(M-|-l)/2 
first-order differential equations (43). In sharp contrast 
to the MA method, the LNA method decouples the com- 
putation of the means from the computation of the co- 
variances. It turns out that the LNA method is substan- 
tially faster than Monte Carlo estimation and can be used 
to provide a rapid assessment of the statistical behavior 
of some Markovian reaction networks (Goutsias, 2006). 
This method has already been used to study biochemical 
reaction networks (Elf and Ehrcnberg, 2003; Hayot and 
Jayaprakash, 2004; Scott et at, 2006; Tao et ai, 2005; 
Tomioka et al., 2004), epidemiological networks (Chen 
and Bokka, 2005), ecological networks (Datta et al., 
2010), social networks (de la Lama et al., 2006), and neu- 
ral networks (Benayoun et al., 2010; Bressloff, 2009). 

F. Macroscopic solution 

For large nonlinear reaction networks, the MA and 
LNA methods can become computationally intractable, 
since evaluation of the covariances requires solving a sys- 
tem of 0{M^) differential equations. If that turns out 
to be the case, then the only option left to character- 
ize the dynamic behavior of the reaction network is in 
terms of DA or population densities by using, for ex- 
ample, the macroscopic (fluctuation-free) system of M 
differential equations given by Eq. (42). As a matter of 
fact, Eq. (41) implies that, for any t > 0, the DA den- 
sity process Zm{t;fl) converges in distribution to Cm{t) 
as f2 — )■ oo. On the other hand, the difference between 
the DA density dynamics predicted by the macroscopic 
system and the MA method grows as ft decreases. In- 
deed, let SC.mit; fi) := fiz{rn; t)/il — C,rnit), where ^z{rn; t) 
is the mean value of the DA process (t) predicted by 
the MA method. Then, from Eqs. (30) and (41) we have 
^-^Wm{t) = 0-V25„(f) _ SCm{t;n) and, since Wm{t) 
is zero mean, we obtain 

6U{t;^) = 4^ E[s™(t)] = o{n-^/^). 



for t > and m G JH. Clearly, for sufficiently large fl, 
6(^rn{t]^) — 0, in which case the macroscopic density 
dynamics obtained by Eq. (42) and the mean density dy- 
namics obtained by the MA method will approximately 
coincide. However, for small fl, 5(rn{t;fl) ^ and 
Eq. (42) may fail to correctly predict the mean density 
dynamics of the DA process. As a matter of fact, it has 
been demonstrated in the literature that, for reaction 
networks with small species populations and appreciable 
stochastic fluctuations (a situation that occurs at small Q 
values), the MA method may reveal behavior that cannot 
be predicted by the macroscopic equation (42) (Gomcz- 
Uribc and Verghese, 2007; Goutsias, 2007; Leonard and 
Rcichl, 1990; McQuarric et al, 1964; Rao and Arkin, 
2003; Srivastava et al, 2002; Thakur et al., 1978; Vel- 
lela and Qian, 2007; Zheng and Ross, 1991). 

Similarly to the DA density process, the population 
density process X{t; Q) converges in distribution, as Q — 

00. to the deterministic process x{t) that satisfies the 
following macroscopic equations: 

^ = E«-^™(XW). (44) 

for i > 0, n e Af, where TTm{x) :— J1^^7rm(f7x), provided 
that these equations are initialized with the same condi- 
tion as the master equation (5). This is clearly true at 
finite times. It is also true in the limit as t ^ oo, pro- 
vided that the macroscopic equations (44) have a unique 
asymptotically stable stationary solution that is indepen- 
dent of the initial state (van Kampen, 2007; Mansour 
et al., 1981). 

G. Remarks 

1. It has been shown by Grima et al. (2011) that, 
for monostable Markovian reaction networks, the LA 
method produces means and covariances that are accu- 
rate to 0(51"^/^) for systems of size fl which are away 
from thermodynamic equilibrium (i.e., systems that do 
not obey detailed balance - see Section VIII) and at least 
accurate to 0{il~^) for systems that are at thermody- 
namic equilibrium (i.e., obey detailed balance). As a con- 
sequence, the LA method will in general result in more 
accurate means and covariances than the LNA method, 
which produces means accurate to 0(17"^/^) and covari- 
ances accurate to 0(f2~^/^). Therefore, Monte Carlo es- 
timation of the means and covariances of the DA and 
population processes based on the LA method will result 
in excellent estimation of the exact values at sufficiently 
large system sizes Q, provided that we use a small time 
step r and a large number L of Monte Carlo samples. 
However, the LA method will produce continuous-valued 
DA trajectories {z{t), t > 0}, which disagrees with the 
fact that these trajectories are integer-valued. 

2. It can be seen from Eqs. (26) and (27) that, in the limit 

as T — > 0"*" , the means and covariances of the approximat- 
ing DA process Z{t) associated with the PA method will 
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approach the means and covariances of the actual DA 
process Z{t) governed by the master equation. Therefore, 
Monte Carlo estimation of the means and covariances of 
the DA process based on the PA method will result in 
excellent estimates of the exact values, provided that we 
use; a sufficiently small time step t and a sufficiently large 
number L of Monte Carlo samples. Note also that, in 
sharp contrast to the Langevin approximation, the Pois- 
son approximation produces DA trajectories {z{t), t > 0} 
that are integer-valued. 

3. When the propensity functions of a Markovian reac- 
tion network satisfy Eq. (39), the macroscopic solution, 
the LNA method, and the MA method provide a hier- 
archy of approximations to the master equation (Ferm 
et al, 2008). At large values of fl, close to the thermo- 
dynamic limit, the macroscopic equations may provide a 
sufficiently accurate description of the reaction network. 
For smaller values of fl, the LNA method will be more 
preferable, whereas, for even smaller values of f2, the MA 
method must be employed. Unfortunately, there is no 
way to determine the range of ft values for which each 
approach is valid. Moreover, for very small values of 
these approximations may not be accurate and Monte 
Carlo simulation methods should be employed instciad. 
It is therefore very difficult to determine a priori which 
method provides accurate estimation of the network dy- 
namics for a given value of O. 

4. For networks governed by the mass-action law, Q usu- 
ally represents the system volume. In this case, the spe- 
cific probability rate constant k of a reaction that involves 
two species will be proportional to with proportion- 
ality factor k; i.e., k = ki}~^ (Gillespie, 1992). This nat- 
urally implies that the frequency of the reaction to occur 
within the infinitcsimally small time interval [t, t + dt) 
will be reduced at a rate that is inversely proportional 
to system volume. As a consequence, the propensity 
function 7r(x; fl) of a reaction Xi + X2 will sat- 
isfy Eq. (39) with f{n) = n, n{x/n) = k{xi/n){x2/n), 

and TT'(x/il) — 0. On the other hand, the propensity 
function of a reaction 2Xi — >• X2 will satisfy Eq. (39) 
with f{n) = n, n{x/n) = {k/2){xi/flf, and 7f'(a;/f2) = 
-{k/2){xi/n). 

5. An important issue associated with the ansatz given 
by Eq. (41) is that, for a given value of fl, the portion of 
the left tail of the Gaussian distribution of the noise com- 
ponent Sm(i) that extends below zero may have appre- 
ciable mass when the standard deviation is large, falsely 
producing negative populations with non-zero probabil- 
ity when such populations are not possible. This problem 
can be taken care of by sufficiently increasing f2, provided 
that the standard deviation of Sm(f) is finite. As a con- 
sequence, justifying the ansatz given by Eq. (41), and 
thus the applicability of the LNA method, requires that 
the covariance matrix Cs{t) is finite. This however may 
not always be true. To see why, note that the unique 
solution of the Lyapunov matrix equation (43) is given 



by (Abou-Kandil et al, 2003) 

€.(*)= f^ait;T)AiT)^lit;T)dT, (45) 
Jo 

where $g (*>''') is an M x M matrix, where M is the 
number of reactions. This matrix satisfies 

^^^^=Gma{t;T), t>T, $^(r;r)=lM, 

with Im being the M x M identity matrix. Since A(t) is a 
diagonal (and thus symmetric) matrix with nonnegative 
elements, Cs{t) will be a symmetric positive semidefi- 
nitc matrix. These properties are required so that Cs{t) 
is a correlation matrix. If the real parts of the eigen- 
values of the Jacobian matrix G(t) are all negative, for 
every i > 0, then, for any fixed r, ^c{t; r) will be finite, 
for every t > t, and limf^oo ^G(i, t) — 0, in which case 
C-s.{t) < 00, for every t > Q!^ This condition is equivalent 
to saying that the solution of Eq. (42) must be asymptot- 
ically stable. As a consequence, lack of asymptotic sta- 
bility of the dynamic evolution of the means of the DA 
density process may result in infinite fluctuations, thus 
violating the ansatz given by Eq. (41) and rendering the 
LNA method invalid. Examples on what can happen in 
this case can be found in van Kampen (1976). 

6. It is interesting to note that A(i) in Eq. (45) is a 
diffusion matrix that tells us how the stochastic proper- 
ties of a reaction network, determined by the propensity 
functions, change at each time point along the mean DA 
trajectory. It turns out that A(t) represents the growth 
of stochastic fluctuations about the mean DA trajectory 
as time progresses. On the other hand, the Jacobian ma- 
trix G(f) produces the dissipation matrix ^dt] t) which 
locally damps the stochastic fluctuations along the mean 
trajectory and squeezes this growth. 

7. The LNA method is not appropriate when the prob- 
ability distributions of the DA and population processes 
are not unimodal, a situation that arises in multistable 
reaction networks (van Kampen, 2007; Qian, 2011). In 
such cases, the approximation can still be applied but 
only during sufficiently short timescales with an initial 
condition that is inside the domain of attraction of an 
equilibrium point of Eq. (42). Moreover, the Gaussian 
nature of the LNA method implies that both processes 
must be continuous- valued. Although this may be ap- 
proximately true for large values of it is not necessarily 
true for smaller values. This is due to the fact that, for 
small fi, the reaction system may contain a small number 
of species interacting through infrequently occurring re- 
actions, in which case Z{t) may take only a small number 
of possible values, at least during an appreciably long ini- 
tial time interval. Note also that, for every value of fl, the 



Note that Kit) is bounded for every t > due to the assumption 
that the propensity functions are analytic. 
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left tail of the Gaussian distribution extends below zero. 
For large fi, the total probability mass over negative DA 
values is negligible and poses no practical problem. Since 
however the covariance matrix of Z{t) is inversely propor- 
tional to n while the mean docs not depend on O, the 
approximation may falsely produce negative DA values 
with appreciable probability when O is small, which may 
erroneously predict negative populations (which is also 
true in the case of the LA and PA methods) , when such 
populations are not possible (Bostani et ai, 2012). 

8. As a consequence of the previous remarks, the LNA 
method provides a reasonable approximation of a Marko- 
vian reaction network exhibiting a single and asymptoti- 
cally stable behavior subject to relatively small stochastic 
fluctuations. By appropriately modifying this method, 
we can also use it to approximate fluctuations around 
mean trajectories ^(t) that produce Jacobian matrices 
G{t) with purely imaginary eigenvalues. In this case how- 
ever the fluctiiations will only be driven by the diffusion 
matrix A(t) and, thus, will grow as a function of time. 
As a consequence, the LNA method can only be used 
over an initial time interval during which the fluctua- 
tions will be sufBciently small so that the ansatz given 
by Eq. (41) can be justified. An appropriate modifica- 
tion of the 0-expansion method can lead in this case to 
approximating fluctuations using the nonlinear Fokker- 
Planck equation or other nonlinear partial differential 
equations (Grima, 2010, 2011; Grima et ai, 2011; van 
Kampen, 1976, 2007). Finally, it may also be possible 
under certain circumstances to modify the LNA method 
to deal with cases in which the eigenvalues of the Ja- 
cobian matrix G{t) have positive real parts [i.e., when 
the mean trajectories ^{t) are unstable [see van Kampen 
(1976, 2007) and Tomita et al. (1974)]. For a recent 
example, see Scott et al. (2006). 

9. A number of moment closure schemes have been pro- 
posed in the literature based on truncating high-order 
central moments or cumulants (Engblom, 2006; Lee et al, 
2009; Matis and Kiffe, 1999, 2002; Nasell, 2003a,b). The 
typical assumption behind these methods is that the so- 
lution of the master equation has negligible high-order 
central moments or cumulants, which can be set equal 
to zero without affecting the mean and covariance dy- 
namics. This assumption is not true in general, with 
the exception of normal random variables whose central 
moments of odd order and cumulants of order > 3 are 
zero. For non-normal random variables, there is an infi- 
nite number of non- vanishing moments or cumulants in 
general (Gardiner, 2010). For example, all cumulants of 
a Poisson random variable are eqiial to the mean value. 
As a consequence, truncation of high-order moments or 
cumulants cannot be easily justified and may lead to a 
non-valid probability distribution (Hanggi and Talkner, 
1980; Hiebeler, 2006; Keeling, 2000). It is worthwhile 
noticing however that low-order cumulants can be used 
to naturally construct univariate and bivariate approxi- 
mations to probability distributions of certain nonlinear 



Markovian processes (Renshaw, 1998, 2000). Moreover, 
it has been suggested by Buice et al. (2010) that an ap- 
propriately defined change of variables that measures the 
deviation of each cumulant from its value under a Pois- 
son assumption (i.e., from the mean) produces a moment 
hierarchy that can be naturally and justifiably truncated 
in the case of neural networks since, in this case, the so- 
lution of the master equation is expected to be near Pois- 
son. Finally, it has been shown by Grima (2012) that, for 
a monostable Markovian reaction network at sufficiently 
large size Q. with at most quadratic propensity functions, 
a normal approximation of the mean of the population 
density process X{t)/Vl (obtained by setting its third- 
and higher-order cumulants equal to zero) is at least as 
accurate as the approximation produced by the macro- 
scopic equations obtained by the Q-expansion method in 
the thermodynamic limit of f2 oo. This approxima- 
tion however may lead to inaccurate covariance dynam- 
ics. On the other hand, a moment approximation scheme 
constructed by setting the fourth- and higher-order cu- 
mulants of the population density process equal to zero 
(thus including third-order central moments in the for- 
mulation) will produce more accurate mean and covari- 
ance dynamics than the normal approximation, provided 
that the system size O is sufficiently large. 

10. The MA method is based on replacing the mean 
value E[a„(Z(t))] by Qm(Mz(i)) + Tm{liz{t)), where 

is a term calculated by solving differential equations 
for higher-order moments [see Eqs. (34) (37)]. When 
the propensity function am(z) is convex, then Jensen's 
inequality implies that Y,[am{Z{t))] > a.m{^[Z{t)]) = 
ctmi/J'zit)). As a consequence, we must have Tm{Hz(t)) > 
0, since E[am{Z{t))] = a„(/Xz(i)) + T^{liz{t)) > 
oim{fji'z{t)). This however may not necessarily be the 
true, in which case the method may result in additional 
errors that can lead to instabilities. To address this 
problem, we can replace E[amiZ{t))] with Q;m(Mz(0) + 
max{0, Tm(/Xz(i))}. In many instances, this simple mod- 
ification results in a more accurate and more stable im- 
plementation of the MA method. 

11. A moment closure scheme used in a particular appli- 
cation must produce moments that satisfy a number of 
necessary conditions. For example, all moments must be 

nonnegative and invariant under permutations. In ad- 
dition, if Zm^{t), Zm^{t), and Zm^it) are mutually un- 
correlated, then we must have E[Zmi{t)Zm2{t)Zm3{t)] 

= ElZ,„,{t)]E[Z„,Jt)]E[Z„,M whereas, if Z,nAt) is 
uncorrelated from Z„i-^{t) and Z^^(t), we must have 

E[ZmAt)ZmAt)Z^sit)] = E[ZmAm[Z^MZmM- 

Moreover, the resulting covariances must define a sym- 
metric positive semi-definite matrix. 

12. The assumptions underlying a given moment closure 
scheme may be inconsistent with the statistics of a partic- 
ular reaction network at hand. In this case, the moment 
equations may have no solution, may produce a num- 
ber of unrealistic solutions, or result in unstable and un- 
bounded solutions (Nasell, 2003a,b). Just like the LNA 
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FIG. 3 The joint probability distributions of the net public and private opinions in the liberal and totalitarian systems of opinion 
formation example of Section IV-I at steady-state computed by the KSA method, (a) The liberal system is characterized by a 
unimodal stationary distribution with its mode located at the point of zero net public and private opinions, (b) The totalitarian 
system is characterized by a bimodal stationary distribution with the modes corresponding to two different totalitarian states; 
one in which a large number of individuals publicly support the ideology of the ruling party, while a small number of individuals 
are privately against this ideology, and one in which a large number of individuals are publicly opposing the current ruling 
ideology, while a small number of individuals privately support it. In the latter case, the ruling party effectively switches its 
current ideology to fit public opinion, thus maintaining itself in power. 



method, normal MA techniques cannot characterize bi- 
modality or highly skewed probability distributions and 
may result in undesirable negative DA and population 
values (Krishnarajah et al., 2005). Finally, much work 
must be done to determine conditions for the stability of 
the differential equations obtained by MA methods. 

H. Example: Opinion formation 

To illustrate and compare the previous methods for 
solving the master equation, we focus on the opinion for- 
mation model discussed in Section III-E. The simplicity 
of this model permits us to solve the underlying mas- 
ter equation using a numerical approach. On the other 
hand, the complexity introduced by the nonlinear nature 
of the propensity functions given by Eq. (9), allows us to 
illustrate some intricate behavior. We consider two pa- 
rameterizations of the model corresponding to a liberal 
democracy and a totalitarian regime (Weidlich, 2006). In 
particular, we assume a social system of 80 individuals 
(in which case L = 40) and set the values of the spe- 
cific probability rate constants associated with the re- 
actions in Eq. (8) to ki = 1/2 day" ^ individual"^ and 
K2 — 1 day ~^ individual"^. 

The main difference between the two social systems 
under consideration is that, in the liberal system, there 
is no pressure inflicted on public opinion with individu- 
als having an affirmative private bias towards the ide- 
ology of the ruling party. On the other hand, there 
is heavy pressure in the totalitarian system inflicted on 
public opinion with individuals having a weakly dissi- 
dent private opinion bias against the ruling ideology. We 
quantify these differences by setting ai = individual"^ 
and as = 1/80 individual"^ in the liberal system, and 



fli = 3/80 individual" , 03 = —1/320 individual" 
in the totalitarian system. In addition, we assume 
that the two systems differ on how strongly privately 
held beliefs infiuence publicly stated opinions and set 
0-2 — 1/80 individual"^ in the liberal system and 02 — 
1/40 individual"^ in the totalitarian system. Finally, we 
consider the condition Xi{0) = ^2(0) = 0, which repre- 
sents completely neutral initial net publicly and privately 
held opinions. 

We can numerically solve the master equation associ- 
ated with the opinion formation model by employing the 
KSA method. This method is more appropriate than the 
IE method, since the DAs of the underlying reactions can 
grow rapidly in this case, whereas the populations Xi{t) 
(net public opinion) and X2{t) (net private opinion) are 
bounded, taking values between —40 and 40. To imple- 
ment the KSA method, we used a Krylov subspace of 
dimension Kq = 40. The joint probability distributions 
of the net public and private opinions in the liberal and 
totalitarian systems at steady-state are depicted in Fig. 3, 
whereas, movies encapsulating the entire dynamic evolu- 
tions of these distributions for a period of 15 days can be 
found in the accompanying Matlab software. Evaluation 
of each solution took about 15 seconds of CPU time on a 
2.20 GHz Intel Core 2 Duo processor running Windows 7. 

In the liberal system, the stationary joint probability 
distribution of public and private opinions is unimodal 
and almost identical to a sampled normal distribution; 
see Fig. 3(a). This distribution characterizes the fact 
that there is no need for an individual to hide her pri- 
vate opinion (thus the high correlation between the pub- 
lic and private opinions). As a consequence, the net 
opinions nearly balance around the origin, which rep- 
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FIG. 4 The stationary marginal probability distributions of 
the net public and private opinions in the liberal system of 
opinion formation example of Section IV-I computed by the 
Krylov subspace approximation method, exact Monte Carlo 
sampling, second-order MaxEnt approximation, and the lineal 
noise approximation method. 

resents an equal number of privately held and publicly 
pronounced opinions for or against the ruling ideology. 
On the other hand, the stationary joint probability dis- 
tribution of public and private opinions in the totalitarian 
system is bimodal; see Fig. 3(b). Here, weakly dissident 
individuals tend to privately disapprove the ruling ide- 
ology, but pressure on public opinion ensures that most 
individuals pubHcly approve this ideology. As a conse- 
quence, dissidence is not strong enough to destabilize 
the totalitarian state and the system will operate around 
the peak located at point (30, —6) with strong net pubhc 
opinion in support of the ruling ideology and a rather 
weak private opinion against this ideology. 

A totalitarian society may move to another peak lo- 
cated at point (—30, 6) with strong net public opinion 
against the ruling ideology and a rather weak private 
opinion in support of this ideology. This situation can be 
easily remedied by the ruling party, which can effectively 
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FIG. 5 The stationary marginal probability distributions of 
the net public and private opinions in the totalitarian sys- 
tem of opinion formation example of Section IV-I, computed 
by the Krylov subspace approximation method, exact Monte 
Carlo sampling, and eighth-order MaxEnt approximation. 

change its ideology to fit the prevailing public opinion, 
thus preventing public upheaval and subsequent removal 
of the totalitarian regime from power. 

Finally, the dynamic evolutions of the joint probability 
distributions provided in the supplemental files demon- 
strate a rapid convergence of the liberal system to its 
stable stationary mode after 3 days and a slower conver- 
gence of the totalitarian system to one of its two modes 
after 14 days. 

Instead of using the KSA method, we can estimate the 
probability distributions of the liberal and totalitarian 
systems by Monte Carlo sampling. To do so, we employ 
4,000 trajectories obtained from the master equation us- 
ing the exact sampling algorithm of Gillespie at a cost 
of about 230 seconds of CPU time. We depict the es- 
timated stationary marginal probability distributions of 
the net public and private opinions in Fig. 4 for the lib- 
eral system, and in Fig. 5 for the totalitarian system. 
For comparison, we also depict the stationary marginal 
probability distributions obtained by the KSA method. 

These results demonstrate the power and weakness of 
Monte Carlo estimation. Monte Carlo sampling is a sim- 
ple and robust approach for solving the master equa- 
tion which, in principle, allows us to compute almost 
any statistical summary of interest with arbitrary preci- 
sion. The probability of a particular event can be esti- 
mated by counting the number of occurrences of the event 
and dividing by the total number of sampled trajectories. 



28 



However, this simple and elegant procedure comes with 
a large computational cost, since the number of trajec- 
tories required to obtain a sufficiently accurate estimate 
is usually very large. Even for estimating the univari- 
ate marginal probability distributions of the net public 
and private opinions, 4,000 iterations does not secim to 
be sufficient, since the symmetric and bimodal nature of 
these distributions can be obscured by estimation errors. 
As a matter of fact, it is often very difficult to know a 
priori how many samples must be used to sufficiently 
estimate the qualitative and quantitative properties of 
a statistical summary of interest or to verify a posteri- 
ori when convergence has occurred. Note however that 
the Monte Carlo method scales far better than the KSA 
method, since trajectories can be sampled from the mas- 
ter equation with a relative ease, even in the case of large 
Markovian reaction networks for which numerical meth- 
ods, such as the KSA or IE methods, cannot be used. 
For large systems, the number of trajectories that can be 
sampled from the master equation in a reasonable time 
will certainly not be adequate for accurately computing 
complex population statistics, such as probability distri- 
butions, but they may be sufficient for estimating certain 
moments (e.g., the means and covariances) of the DA and 
population processes with a desired precision. 

To ameliorate the computational burden of exact 
Monte Carlo sampling, we can employ a Langevin or 
Poisson approximation. Sampling 4,000 trajectories from 
the master equation governing the opinion formation 
model using the LA method took 30 seconds of CPU 
time, whereas, drawing the same number of samples us- 
ing the PA method required 75 seconds of CPU time. The 
estimation results are very similar to the ones obtained 
by exact Monte Carlo sampling (data not shown). The 
LA method however does not retain the integer-valued 
nature of the net opinion trajectories, which can be an 
issue when accurate simulation of these trajectories is 
desired. 

By using the PA method, we can draw 4,000 sam- 
ples from the master equation and use these samples to 
estimate the first K moments of the opinion trajecto- 
ries under consideration. We can then use the MaxEnt 
method discussed in Section IV-C-5 to derive an analyt- 
ical approximation of the marginal probability distribu- 
tions of the net public and private opinions. The resulting 
sampled probability density functions depicted in Figs. 4 
and 5 clearly show the potential of the MaxEnt method 
for correctly estimating the stationary marginal distri- 
butions in the opinion formation model. It took about 
75 seconds of CPU time to compute these results, which 
have been obtained by using the Matlab code developed 
by Mohammad- Djafari (1991). For the democratic sys- 
tem, we only need to estimate the first- and second-order 
moments by Monte Carlo. This is due to the fact that 
the true stationary distributions are almost sampled nor- 
mal. On the other hand, to obtain sufficiently accurate 
approximations of the stationary marginal distributions 
in the totalitarian system, we need to estimate the first 



eight moments by Monte Carlo. The main advantage of 
MaxEnt is its ability to provide a relatively accurate ap- 
proximation of marginal probability distributions by us- 
ing appreciably fewer sample trajectories than the ones 
required by Monte Carlo in order to achieve a similar 
level of estimation accuracy. 

Application of the LNA method for solving the mas- 
ter equation associated with the totalitarian system is 
not possible due to the bimodal nature of the stationary 
joint probability distribution (see Remark IV-G-7). This 
method however can be used in the case of the liberal 
system by choosing the system size parameter il to be 
equal to the "size" L of the net public or private opin- 
ions. Evaluation of the solution obtained by the LNA 
method took only 0.35 seconds of CPU time. The result- 
ing Gaussian probability density function approximates 
well the stationary solution found by the KSA method. 
The sampled computed stationary marginal probability 
distributions of the net public and private opinions arc 
depicted in Fig. 4. If there were more individuals in 
the system (i.e., for larger values of ft), then the LNA 
method could produce a more accurate result. On the 
other hand, a significantly smaller number of individu- 
als may dramatically reduce the accuracy of the method, 
since the statistical properties of the system may appre- 
ciably deviate from normality. Despite its clear compu- 
tational advantage, use of the LNA method is hampered 
by the absence of a strategy to effectively determine for 
which values of fl the resulting normal approximation is 
accurate. 

Finally, use of the MA method is much easier for the 

liberal system than the totalitarian system, since the to- 
talitarian system requires at least eighth-order moments 
to sufficiently characterize its bimodal stationary distri- 
bution. For simplicity, we will therefore focus on the 
liberal system. Note that, due to the exponential nature 
of the propensity functions, given by Eq. (9), their effect 
on the moment equations can persist through infinitely 
many derivatives, which can make the task of finding an 
appropriate moment closure scheme very difficult. How- 
ever, since the stationary joint probability distributions 
of the net public and private opinions approximate well 
a sampled Gaussian distribution (see Fig. 4), we may be 
able to approximate the means and covariances by using 
the normal MA scheme given by Eqs. (36) and (37). In 
Fig. 6, we depict the means (solid red lines) and the ±1 
standard deviations (dashed lines) of the net public and 
private opinions, obtained with the normal MA method. 
Implementation of this method took a mere 0.6 seconds 
of CPU time. For comparison, we also depict the cor- 
responding moments and standard deviations obtained 
with the KSA method. Clearly, the normal MA method 
produces unsatisfactory results. 

The main culprit here is the fact that Eqs. (36) 
and (37) were derived for quadratic propensity func- 
tions whose higher-order derivatives vanish, which along 
with the Gaussian assumption results in a decoupling 
of the means and covariances from higher-order central 
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FIG. 6 The means (solid red lines) and the ±1 standard de- 
viations (dashed lines) of the net public and private opinions 
in the liberal system of opinion formation example of Sec- 
tion IV-I, obtained by the Krylov subspace and normal mo- 
ment approximation methods (with and without correction 
using Jensen's inequality). 

moments. In the democratic system, the Gaussianity 
assumption is approximately valid but errors accumu- 
late, since Eqs. (36) and (37) neglect to account for non- 
vanishing higher-order derivatives of the propensity func- 
tions. We can however improve the closure scheme by 
using Jensen's inequality (see Remark IV-G-10), due to 
the convexity of the propensity functions. The results 
depicted in Fig. 6 clearly demonstrate the effectiveness 
of this correction. 

V. MULTISCALE METHODS 

As we discussed earlier in this review (see Section IV- 
C-6), the reactions in a Markovian reaction network may 
occur at different time scales, with slow reactions oc- 
curring infrequently and fast reactions firing numerous 
times between successive occurrences of slow reactions. 
This may appreciably increase the computational effort 
required to sample the master equation by Monte Carlo, 
which can make analysis of Markovian reaction networks 
very difhcult to perform in practice. In this section, we 
discuss methods available to address this problem. The 
main idea is to eliminate the fast reactions by approx- 
imating the master equation with one that consists of 
only slow reactions. 



A. Partitioning approximation 

In many cases of interest, it is not important to know 
the detailed activities of fast reactions, since the dynamic 
evolution of the state of a Markovian reaction network 
may be mostly determined by the slow reactions. If that 
is true, we may be able to approximate the master equa- 
tion by one that consists only of slow reactions. If a 
sufficiently accurate approximation of the master equa- 
tion can be found in terms of slow reactions, then it can 
be used to appreciably reduce the computational com- 
plexity associated with Monte Carlo sampling. This is 
due to the fact that sampling slow reactions is appre- 
ciably more efficient than sampling fast reactions. This 
idea has led to the development of techniques for elimi- 
nating fast reactions, known as multiscale or partitioning 
approximation methods (Ball et al., 2006; Bostani et al., 
2012; Cao et ai, 2005b; Chevalier and El-Samad, 2009; 
Cotter et al, 2011; Crudu et al, 2009; E et al, 2005, 
2007; Gomez-Uribe et ai, 2008; Goutsias, 2005; Hasel- 
tine and Rawhngs, 2002, 2005; Kang and Kurtz, 2012; 
Mastny et ai, 2007; Pahlajani et ai, 2011; Puchalka and 
Kierzek, 2004; Rao and Arkin, 2003; Sails and Kaznessis, 
2005; Samant and Vlachos, 2005; Thomas et ai, 2012; 
Zeron and Santillan, 2010). 

To illustrate the main steps underlying most multi- 
scale approximation schemes, let us assume that the 
first Ms reactions of a Markovian reaction network are 
slow, whereas, the remaining Alf — M — Ms reactions are 
fast. We can then decompose the DA process Z(t) into 
two components, Zs{t) and Zf(t), with the first compo- 
nent corresponding to the slow reactions and the second 
corresponding to the fast reactions. In this case, Eq. (4) 
leads to the following master equation (Goutsias, 2005): 



{zs;t)pz^{zs;t)y (46) 



dt 



me Ms 



for t > 0, with 



2/ 



{Zs,Zf)p^^l^^{Zs,Zf;t), 



(47) 



for m G Ais , where pz^ (Zg ; t) is the marginal proba- 
bility distribution over the DAs of the slow reactions, 
Pzf\zsizs,Zf]t) is the conditional probability of the fast 
DAs at time t given the DAs of the slow reactions, 
Ms ■= {1,2,..., Ms}, and is a vector comprised of 
the first Ms rows of e,„ (the m-th column of the M x M 
identity matrix). 

According to Eq. (46), the DAs of the slow reactions 
follow a master equation that is similar to the one gov- 
erning the entire Markovian reaction network, albeit with 
time-varying propensities. The fast reactions exert their 
influence on the slow reactions through the propensity 
functions given by Eq. (47), which are computed as the 
conditional means of the original propensity functions 
given the DAs of the slow reactions. As a consequence, if 
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we can evaluate the mean propensity functions am {zs ; t) , 
m G A^s, then we can efficiently simulate the stochas- 
tic evolutions of the DAs of the slow reactions by us- 
ing the exact Gillespie algorithm, or any other appropri- 
ate technique, modified to account for the fact that the 
propensity functions are now time-dependent (Haseltine 
and Rawlings, 2002; Rao and Arkin, 2003). Moreover, 
given that Zs{t) = Zg, we can estimate the population 
process Xn{t) by using the optimum mean-square esti- 
mate x„(t), given by [recall Eq. (3)]: 

Xn{t) = '&[Xn{t)\Zs{t)=Zs\ 

= Xo,n + XI SnmZm{t) +^ Snml^z{m;t, Zg), (48) 
m^Ms m&Mf 

for n e J\f, where fj,z{m;t, Zg) := E[Z„i(i) | Zs{t) = z^], 
for m€ Mf := {Mg + l,Ms + 2,..., M}, is the mean DA 
of the m-th fast reaction at time t, given the state Zs{t) 
of the slow reactions at t. Hence, we can replace the 
Markovian reaction network by a "slow" Markovian re- 
action subnetwork whose DA process is governed by 
Eqs. (46), (47) and whose population process is governed 
by Eq. (48). 

Calculating the propensity functions (^s ; ^) , "ti € 
A^s, and the means fizimit, z^) , m € Aif , requires knowl- 
edge of the conditional probability Pj, j, (Zg , Z/ ; t) . It can 
be shown that, within the coarse time scale, the dynamic 
evolution of this probability is approximately governed 
by the following master equation (Goutsias, 2005): 

dPzf\zA'^s,Zf;t) 
dt 

= ^ [am{Zs,Zf - e^)Pzj {Zs ,Zf-e^;t) 

m&Mf 

-am{zs,Zf)pzf\zS^s^^f^t)}^ (49) 

for t > 0, where e„ is a vector comprised of the last Aif 
rows of Cm, which includes only the fast reactions. This 
equation is derived by assuming that, within successive 
firings of slow reactions, the "slow" Markovian reaction 
subnetwork behaves like one with propensity functions 
that are not appreciably larger than zero, since it is very 
unlikely that a slow reaction will occur within that time 
scale. 

It turns out that most multiscale approximation 
schemes currently available in the literature follow a 
similar approach by decomposing a Markovian reaction 
network into "slow" and "fast" Markovian subnetworks 
based on appropriately chosen variables. Although the 
method discussed above uses the DAs of individual reac- 
tions (Goutsias, 2005; Haseltine and Rawlings, 2002; Salis 
and Kaznessis, 2005), similar methods can be constructed 
using the net DAs of reversible^'^ reactions (Haseltine 



A reaction is called reversible if it can occur in both directions, 



and Rawlings, 2005), or the populations of the underly- 
ing species (Cao et ai, 2005b; Chevalier and El-Samad, 
2009; Gomez-Uribe et al, 2008; Rao and Arkin, 2003; 
Samant and Vlachos, 2005; Zcron and Santillan, 2010). 

Unfortunately, solving the master equation of the 
"fast" reaction subnetwork [e.g., Eq. (49)] is as difficult as 
solving the master equation of the entire network. More- 
over, evaluating the propensity functions of the "slow" 
subnetwork requires Monte Carlo estimation in general, 
which adds to computational complexity. To address 
these issues, a number of different approaches have been 
proposed in the litcrattire. ba.scxl on the techniques dis- 
cussed in Section IV. For example, it has been assumed 
that, within successive firings of slow reactions, the fast 
reactions rapidly reach a stationary state whose proba- 
bility does not depend on time t (Cao et al, 2005b; E 
et al, 2005, 2007; Gomez-Uribe et al, 2008; Haseltine 
and Rawlings, 2005; Rao and Arkin, 2003; Samant and 
Vlachos, 2005; Zeron and Santillan, 2010). In this case, 
we can set the right-hand-side of the master equation of 
the "fast" reaction subnetwork equal to zero and use a nu- 
merical technique (see Sec. IV-B) to calculate the desired 
stationary conditional probability of the "fast" variables 
given the "slow" variables. We can then evaluate the 
propensity functions of the "slow" reaction subnetwork 
either by direct summation (Cao et al, 2005b; Haseltine 
and Rawlings, 2005), if computationally feasible, or by 
Monte Carlo estimation (Haseltine and Rawlings, 2005). 

Numerically solving the master equation of the "fast" 
reaction subnetwork may not be easy, especially for large 
subnetworks. Moreover, evaluating expectations by di- 
rect summation or Monte Carlo estimation can be com- 
putationally demanding. The main difficulty however 
with the previous approach is to verify that the "fast" re- 
action subnetwork reaches a stationary state, since there 
might be only a short induction time between succes- 
sive firings of slow reactions during which convergence to 
steady-state may not occur. 

We can avoid the previous problems by sampling the 
master equation of the "fast" reaction subnetwork using 
the exact Gillespie algorithm and evaluate the propensity 
functions of the "slow" reaction subnetwork by Monte 
Carlo estimation (E et al, 2005, 2007; Samant and Vla- 
chos, 2005). Although this strategy is very general, re- 
quiring no additional assumptions other than the ones 
leading to Eqs. (46), (47), and (49), the embedded esti- 
mation step may require a large number of Monte Carlo 
samples, which can substantially increase computations. 

Exact sampling of the master equation of the 
"fast" reaction subnetwork can be replaced by the 
Langevin (Haseltine and Rawlings, 2002, 2005; Salis and 
Kaznessis, 2005) or the Poisson (Puchalka and Kierzek, 
2004) approximations. Although this can speed-up sam- 



arbitrarily labeled as "forward" and "backward," with nonzero 
probability; otherwise, the reaction is called irreversible. 
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pling of the "fast" variables, it will not eliminate the need 
for evaluating the propensity functions of the "slow" re- 
action subnetwork using Monte Carlo estimation. Due 
to its discrete nature, the Poisson approximation may be 
more preferable than the Langevin approximation. How- 
ever, for both approximations to be valid and computa- 
tionally efficient, it is necessary to determine a value for 
the leaping parameter r that is as largo as possible while 
still ensuring that occurrences of fast reactions within 
a time interval [t, t + t) do not appreciably affect the 
propensity functions. Intuitively speaking, finding such 
value may be possible due to the assumed futility of the 
fast reactions. In practice however this may not be easy. 
Note finally that the expected number of occurrences of 
each fast reaction during [t, t + t) will be much larger 
than one, a condition that is required for the Langevin 
approximation to be valid. 

The propensity ftmctions of the "slow" reaction sub- 
network can be approximated by using Eq. (47) and a 
Taylor series expansion, such as the one given by Eq. (33) , 
of the propensity functions of the entire network around 
the means of the "fast" variables. This can be done by 
evaluating the conditional moments of the "fast" vari- 
ables, given the values of the "slow" variables, thus elim- 
inating the need for Monte Carlo estimation. If the 
propensity functions of the "slow" reaction subnetwork 
depend only linearly on the "fast" variables, then we only 
need to calculate the conditional means of these variables. 
On the other hand, if the "slow" propensity functions de- 
pend quadratically on "fast" variables, then we also need 
to calculate the conditional covaxiances. We can per- 
form these calculations by employing a moment approx- 
imation scheme applied on the master equation of the 
"fast" reaction subnetwork (Cao et al, 2005b; Chevalier 
and El-Samad, 2009; Gomez-Uribe et al., 2008; Goutsias, 
2005; Haseltine and Rawlings, 2002, 2005; Rao and Arkin, 
2003). The accuracy of this approach however depends 
on the degree of nonlinearity of the network propensity 
functions in terms of the "fast" variables and the partic- 
ular moment approximation scheme used. 

Most multiscale approximation methods developed so 
far are based on a clear separation between fast and 
slow reactions. In reality however this may not be pos- 
sible. For this reason, it may be more appropriate to 
develop techniques that involve more than two separate 
time scales. We refer the reader to E et al. (2005, 2007) 
and Harris and Clancy (2006) for two promising tech- 
niques along this direction. 

We should finally point out that a few alternative mul- 
tiscale approximation schemes have been proposed in the 
literature, namely two techniques related to the finite 
state projection method (Peles et al., 2006; Pigolotti and 
Vulpiani, 2008), a method based on separating species in 
terms of their variance (Hellander and Lotstcdt, 2007), a 
technique based on an adiabatic approximation using a 
stochastic path integral (Sinitsyn et al., 2009), and a rig- 
orous and versatile technique based primarily on stochas- 
tic equations determining the dynamic evolution of the 



population process itself (Ball et al., 2006; Kang and 
Kurtz, 2012). Although promising, those methods have 
only been applied to very small reaction networks. At 
this point, it is not clear how they will perform when 
dealing with larger and more realistic networks, nor have 
they been sufficiently compared to the other approaches 
discussed in this section. 

B. Example: Transcription regulation 

To illustrate the effectiveness of a multiscale approxi- 
mation method for solving the master equation, wo con- 
sider here a biochemical reaction network comprised of 
six molecular species that interact through the following 
ten reactions: 



reaction 1 




Xi 


XI + X2 


reaction 2 




2X2 


^ X3 


reaction 3 




X3 


2X2 


reaction 4 




+ Xi 


^ X5 


reaction 5 




X5 


^ X3 + X4 


reaction 6 


X3 


+ X5 


^ Xe 


reaction 7 




Xe 


^ X3 + X5 


reaction 8 




x^ 


^ X1 + X5 


reaction 9 




Xi 


^ 


reaction 10 : 


X2 


^ 0. 



This reaction network was originally proposed by Gout- 
sias (2005) and can serve as a model for a particular type 
of transcription regulation in single cells. As a matter 
of fact, reaction 1 can be used to model the translation 
of an mRNA molecule Xi into a protein molecule X2, 
whereas reactions 2 and 3 can be used to model dimer- 
ization of X2 into X3. On the other hand, reactions 4- 
7 can be used to model the binding of dimer X^ on a 
gene X4, assuming that the promoter of this gene has 
two binding sites for X3, whereas, reaction 8 can bo used 
to model the transcription of X4 to mRNA molecules Xi , 
assuming that this process occurs when the promoter 
of X4 is only bound by one dimer X3. Finally, reac- 
tions 9 and 10 model degradation of the mRNA and 
protein molecules Xi and X2, respectively. Here, we 
slightly simplify the original model by assuming that the 
cell's volume remains fixed at 10^^^ liters (1). In this 
case, the specific probability rate constants are set to 
Ki = 0.043 s-i, K2 = 0.083 moles • -s-i, K3 = 0.5 s^^, 
K4 = 0.0199 moles • • s-\ K5 = 0.4791 s-\ Kg = 
1.9926 X 10"'' moles • ■ s^S K7 = 8.7658 x lO"^^ y-i^ 
Kg = 0.0715 s-\ Kg = 0.0039 s-\ and mo = 0.0007 s'^, 
in agreement with the values used by Goutsias (2005). 
Finally, we initiahze the system by setting Xi{0) = 0, 
X2{0) = 2, X3(0) = 4, ^4(0) = 2, and ^5(0) = ^6(0) = 
0; i.e., we assume that the system contains initially two 
mRNA molecules, two copies of the same gene, and four 
dimer s. 

Despite the modest size of the previous network, sim- 
ulation using exact Monte Carlo sampling of the mas- 
ter equation is computationally intensive. It took about 
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2 hours and 15 minutes of CPU time on a 2.20 GHz In- 
tel Core 2 Duo processor running Windows 7 to obtain 
2,000 samples of the population dynamics during a pe- 
riod of 35 minutes. This serious inefficiency is due to the 
stiffness of the system caused by the reversible reactions 
associated with the dimerization of protein X2 (i.e., reac- 
tions 2 and 3) being much faster than the remaining reac- 
tions. As a consequence, exact Monte Carlo sampling is 
forced to spend a substantial amount of time simulating 
the occurrences of these two reactions. Unfortunately, we 
cannot appreciably reduce computational effort by using 
the Poisson approximation since, for accurately solving 
the master equation, stiffness constrains the leaping pa- 
rameter r to take a very small value thus deeming this 
approximation method computationally comparable to 
exact sampling. Dimerization however is reversible and 
occurs on a much faster timescale than the other reac- 
tions. As a consequence, we expect its effect to largely 
cancel out. Therefore, faithful simulation of dimerization 
may not be necessary. 

If we set Ms = {1, 4, 5, 6, 7, 8, 9, 10} and 7W/ = {2, 3}, 
then the "slow" subsystem, comprised of the reactions 
iTL AAs, will be characterized by the master equation (46) 
with propensity functions given by [recall Eq. (47)] 

ai{zs]t) = Ki(z8 - Zg), 

ai(zs;t) = K4[4 - Z4 + Z5 - Z6 + Z7 + Mz(2; t,Zs) 
-^lz{'i■,t,Zs)]{2- Zi + Z5), 

a^{Zs]t) = K5(Z4 - Z5 - Z6 + Z7), 

aQ{Zs]t) = Z4^ + zq + zt + fiz{2;t,Zs) 

-Hz{i;t,Zs)]{z4 - zz~ ze + Z7), 
a7{zs;t) = Kr{ze - zj), 
as{zs; t) = Kfi{zi - z^- z(i + Z7), 
a9(zs;i) = K9(z8 - zg), 

aio {Zs ; i) = Kio [2 + zi - zio - 2/iz (2; i , z^) + 2/1^ (3; t , z^)] , 

where /Xz(2;t,Zs) and /iz(3;i,Zs) are the mean DAs of 
the two fast reactions 2 and 3, respectively. Therefore, 
to calculate these propensities, we need to compute the 
difference /Xz(2; t, Zg) — jizi^] t, Zg). By assuming that the 
"fast" reaction subsystem of the two dimerization reac- 
tions rapidly reaches equilibrium within successive occur- 
rences of slow reactions, we can show that [see Goutsias 
(2005)] 



^z(2;t,Zs) - ^izi^;t,Zs 



where 
A{zs) 



(50) 



1.5 + zi - zio + (K3/2K2) 



B(z,) = 0.25(1 + zi - zio)(2 + zi - 
- (k3/2k2)(4 - Z4 + Z5 



zw) 
Z6 ^ 



(51) 



Z7 



As a consequence, we can solve the master equation (46) 
of the "slow" reaction subsystem without having to 
solve the conditional master equation (49) of the "fast" 
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FIG. 7 Means (solid lines) and ±1 standard deviations 
(dashed lines) of the population processes in the transcrip- 
tion regulation network example of Section V-B obtained by 
exact Monte Carlo sampling (blue lines) and multiscale ap- 
proximation (red lines). 

subsystem, and estimate the population process by us- 
ing Eq. (48), which depends only on the difference 
Aiz(2;t,z,)-Aiz(3;t,z,), and Eqs. (50), (51). 

It took less than a minute (52 seconds) of CPU time to 
draw 2,000 Monte Carlo samples from the master equa- 
tion of the "slow" reaction subsystem (as compared to 
135 minutes of CPU time required by exact Monte Carlo 
sampling). The mean and ±1 standard deviation dy- 
namics of the underlying population processes obtained 
by the two methods are depicted in Fig. 7. These re- 
sults clearly show that an appropriately derived multi- 
scale approximation of a stiff Markovian reaction network 
can lead to dramatic improvements in computational ef- 
ficiency while producing a relatively accurate approxi- 
mation of the population dynamics. It turns out that 
the relatively large transient errors in the population dy- 
namics of X4 and depicted in Fig. 7 are due to in- 
correctly computing the net DA 2'(4) — 2'(5) of reac- 
tions 4 and 5 (binding and unbinding of dimer X-^ on 
the promoter of gene ^4), which is a consequence of the 
imposed adiabatic approximation of the DAs of the fast 
reactions 2 and 3 (dimerization) through their mean val- 
ues. Although this approximation also affects the accu- 
racy of the remaining population dynamics, the resulting 
approximate dynamics track the ones computed by exact 
Monte Carlo sampling sufficiently well. 
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VI. MESOSCOPIC (PROBABILISTIC) BEHAVIOR 

An important goal when studying Markovian reaction 
networks is to investigate the existence, uniciueness, and 
stability of a stationary solution of the underlying mas- 
ter equation and derive mathematical properties of the 
dynamic behavior of the probability distribution of the 
system state. This can be done by using a mesoscopic de- 
scription of the network in terms of the population prob- 
abilities {px{x;t),x e X}, for t > 0. To avoid math- 
ematical subtleties, which are outside the scope of this 
review, we assume here that the cardinality of the pop- 
ulation state-space X is finite. Most results however can 
be extended to the case of countable state-spaces. 

To derive a stationary solution of the master equa- 
tion (5), we must solve the system of K linear equations 
Pp = 0; recall Eq. (20). Since the elements of each col- 
umn of matrix P add to zero, its rows are linearly depen- 
dent and, therefore, the rank of P will be less than K. As 
a consequence, Pp = will have at least one nontrivial 
solution. Unfortunately, this result does not tell us how 
many nontrivial solutions exist and which ones are valid 
probability distributions; i.e., which solutions satisfy the 
necessary constraints 



< P/c < 1, for fc = 1, 2, . . . , i^, and 



K 

E 

fe=i 



Pk = 1- 



In the following, we first consider the dynamic behav- 
ior of an irreducible Markovian reaction network. This 
type of network is defined by the property that, for any 
pair {x,x') of population states, there exists at least one 
sequence of reactions that takes the system from state x 
to state x' - these states are said to be communicating. 
By using a simple graph theoretic analysis and Kirch- 
hoff's theorem, Schnakenberg (1976) has shown that an 
irreducible Markovian reaction network converges to a 
unique probability distribution p at steady-state, which 
does not depend on the initial probability distribution 
p(0), such that < p < 1^^ [see also van Kampen (2007)]. 
As a consequence, in an irreducible Markovian reaction 
network, the population process can take any value in X 
at steady-state with nonzero probability. 

On the other hand, the theory of systems of ordi- 
nary differential equations with constant coefficients im- 
plies that, for a given initial probability distribution p(0), 
Eq. (20) is satisfied by a unique probability distribution 
p{t), t > 0, which is analytic for all < t < oo. Since the 
elements of each column of matrix P add to zero, 



d[e^p(t)] T dp{t) 



dt 



dt 



: e^Pp(i) = 0, 



where e is a ii' x 1 vector with all its elements being equal 
to one. This result, together with the fact that e-^p(O) = 



a < V < b denotes that each element of vector v satisfies this 
relation. 



1, implies e^p{t) = 1, for all t > Q. Unfortunately, it is 
not clear whether < p{t) < 1, for every t > 0. It turns 
out however that, for an irreducible Markovian reaction 
network, < p{t) < 1, for every t > (Schnakenberg, 
1976). 

If Afe, k = 1,2, . . . ,K, are the eigenvalues of matrix P, 
with corresponding right and left eigenvectors rk,lk, k = 
1,2, ... ,K, respectively, then the solution to Eq. (20) is 
given by (Moler and van Loan, 2003) 



K 



p{t) = exp (Pt)p(O) = ^ Cfe rfe e^-* 



(52) 



fe=i 



for < t < DO, where we assume here that the eigenval- 
ues of P have the same algebraic and geometric multiplic- 
ity, an assumption satisfied by many Markovian reaction 
networks. In this case, the right and left eigenvectors are 
biorthogonal (i.e., Ij^rk' = 0, for every k ^ k'), which im- 
plies that the constants are given by = l^p{0)/l^rk. 
As a consequence, we can use the eigenvalues and eigen- 
vectors of P to analytically specify the entire mesoscopic 
behavior of a Markovian reaction network. See Keeling 
and Ross (2008) and Vellela and Qian (2009) for appli- 
cation of Eq. (52) to problems in epidemiology and com- 
putational biochemistry. 

Note that Eq. (52) and the fact that a non-trivial sta- 
tionary solution always exists imply that at least one 
eigenvalue of P must be zero. For an irreducible Marko- 
vian reaction network, matrix P has only one zero eigen- 
value, with the remaining K — 1 eigenvalues having neg- 
ative real parts (Schnakenberg, 1976). If we therefore as- 
sume that Ai = 0, then Eq. (52) implies that the station- 
ary distribution will be given by p = ri/||ri||, where ri 
is the eigenvector corresponding to the zero eigenvalue 
and ||r|| is the £i-norm of vector r. It turns out that the 
solution p(i), t > 0, ol Eq. (20) is asymptotically stable 
with respect to p, in the sense that 



lim D\p{t),p] = 0, 

t—^oo 



where 



K 



D\p,q] := ^Pfcln 



Pk 

Qk 



(53) 



is the KuUback-Leibler distance between the two prob- 
ability distributions p = {pk,k = 1,2,...,K} and 
q = {qk,k = 1,2,...,/!'}. As a matter of fact, 
dD\p{t),p]/dt < 0, where equality is archived only at 
steady-state. 

To summarize, for a given initial probability vector 
p(0), the master equation associated with an irreducible 
Markovian reaction network has a unique and strictly 
positive solution < p{t) < 1, < t < oo. This solution 
is analytic for all < t < oo, converges to a stationary 
distribution < p < 1 that does not depend on the initial 
probability distribution p(0), and is asymptotically stable 
with respect to p. 
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It is not in general easy to check whether a Marko- 
vian reaction network is irreducible. However, we often 
assume that a given Markovian reaction network is com- 
prised of only reversible reactions (reactions which can 
occur in both directions with nonzero probability). This 
is a plausible assumption since, in principle, a transi- 
tion between two physical states can occur in the reverse 
direction as well. In this case, and after appropriately 
ordering the states, matrix P can be cast into a block di- 
agonal form with diagonal elements p(^),p(^\ . . . ,p(''), 
for some J, where each submatrix P^-'^ is irreducible 
(when J = 1, matrix P is itself irreducible). The re- 
sulting Markovian reaction network is said to be com- 
pletely reducible (van Kampen, 2007). In this case, the 
original Markovian reaction network can be decomposed 
into J non-interacting subnetworks with non-overlapping 
state-spaces, which can be treated independently of each 
other. Each reaction subnetwork is characterized by 
unique dynamic and stationary solutions p^^^'lt), p^^\ 
j ~ 1, 2, . . . , J, that satisfy the aforementioned proper- 
ties. However, the dynamic and stationary solutions of 
the original master equation are determined by the ini- 
tial condition at time t = 0. If the master equation is 
initialized with a population vector in the state-space of 
the j-th subnetwork, then its dynamic and stationary so- 
lution will be given by 



p(-'')(t) 



and 



,0) 



respectively, where p'^^\t) depends on the initial condi- 
tion and p'^' does not. 

A question that arises at this point is what hap- 
pens when the Markovian reaction network contains irre- 
versible reactions and matrix P is not irreducible. To get 
an idea, let us assume that, after appropriately ordering 
the states, 



p(i) jW 
T(2) 



where P^^) and T^^) are square matrices, P^-"^^ is irre- 
ducible, and at least one element of each column of T^^^ is 
strictly positive. Note that the nonzero elements of T^^^ 
correspond to nonreversible reactions. The associated 
Markovian reaction network is said to be incompletely re- 
ducible (van Kampen, 2007). If we denote byp^^^(t) and 
p^'^\t) the probability distributions of the state vectors 
at time t, determined by the partition of the state-space 
suggested by matrix P, then the master equation results 



in the following two differential equations: 

*^=T(V^)(i). 

Clearly, one can solve the second equation independently 
from the first to obtain 

p(2nt)=exp{T(2)f}p(2)(0). 

On the other hand, the dynamic behavior of p^^^ is 
now driven by p(^^(t), unless p^'^\0) = 0, in which case 
p(i)(t) = exp{p(i)i}p(i)(0). Note however that 



dt 



dt 



= _eTTWpW(f) <0, 

provided that p'^' (t) ^ 0, for every t > 0, since the el- 
ements of each column of matrix P add to zero and we 
have assumed that each column of matrix T^^^ contains 
at least one element that is strictly positive. Therefore, 
p^^\t) asymptotically becomes zero as t ^ oo. As a 
matter of fact, p^^\t) assigns probability mass over the 
transient states of the Markovian reaction network, as 
opposed to p'^^ (t) that assigns probabiUty mass over the 
persistent states. In this case, and when matrix P^^^ is ir- 
reducible, the stationary solution of the master equation 
governing an incompletely reducible Markovian reaction 
network will be unique and given by the probability vec- 
tor 



P 



(1) 



where p(i) is the (unique) solution of the linear system 
of equations pf-'^'p = 0. 

In general, the population states in a Markovian reac- 
tion network can be classified into two distinct groups: 
transient and persistent. These states can be uniquely 
partitioned into non-overlapping sets T and Pj, j = 
1, 2, . . . , J, where T contains all transient states and Pj, 
j = 1, 2, . . . , J, are irreducible sets containing persistent 
states with the additional property that, for every j ^ j' , 
each state in Pj does not communicate with any state 
in Pj/ . By appropriately ordering the states, we can write 
matrix P in the form 



p(i) 

p(2) 



t(i) "1 
T(2) 



• • • Pf"') T^"') 
• • • T 
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where P^^^ is a square irreducible matrix that charac- 
terizes how probabiUty mass is dynamicahy distributed 
among the persistent states in Pj, T'^-'-' is a matrix that 
tells us how probability mass is transferred from the tran- 
sient states in T to persistent states in Pj, and T is a 
square matrix that characterizes how probability mass is 
dynamically distributed among the transient states in T. 
In this case, if the Markovian reaction network is ini- 
tialized by a persistent state in Pj, then the stationary 
solution will be given by the probability vector 





93 



pU) 







where p^-'^ is the unique stationary distribution of the j-th 

irreducible Markovian reaction subnetwork characterized 
by matrix P^j\ However, if the network is initialized 
with the i-th transient state in T, then the stationary 
distribution (which now depends on i) will be given by 
a convex combination of the stationary distributions gj 
above, with mixing coefficients i.e., we have that 

,7 

Pi = ^1^139 h (54) 

where 

J 

Hij > and ^^Mij — 1- 

i=i 

As a matter of fact, Eq. (54) simply expresses the fact 
that the probability of a Markovian reaction network ini- 
tialized with the i-th transient state in T to reach a per- 
sistent population state x in Pj at steady-state equals 
the probability fiij that the system will reach a persistent 
state in Pj at steady-state multiplied by the probability 
that this state will be x. It can be shown that 

i'eTj'ePj 

where [T'^-'^Jj/i/ is the clement of matrix T^^^ and 

[T~^]j'j is the {i' , i) element of the inverse of matrix T. 

To summarize, a fundamental property of the master 
equation (5) associated with a Markovian reaction net- 
work is that, when this equation is initialized with a per- 
sistent state, its solution converges to a unique stationary 
distribution that assigns positive probability only to the 
persistent states that communicate with the initial state. 
On the other hand, if the Markovian reaction network 
is initialized with a transient state, then its stationary 
distribution will be a convex combination of the distinct 
stationary distributions obtained by initializing the sys- 
tem with persistent states chosen from each individual 
irreducible set. 



VII. POTENTIAL ENERGY LANDSCAPE 

To better understand what might happen at steady- 
state, let us assume that the master equation (5) has 
a unique stationary solution Py^{x) :~ limt^ao Px{x',t) 
that is independent of the initial state. The probabil- 
ity distribution Px{x; t) of the population density process 
X{t] n) = X{t)/n will be given by ps(x; t) = Vlp^inx; t) 
and will depend on the size parameter Q. in general. Let 
us define the function 



V{x; n) := 



1 1^ Px{^) 

o Px{x^y 



(55) 



where [x) :~ lim(_>.oo Px (5; i) is the steady-state distri- 
bution of the population density process and 5* is a state 
at which the stationary probability distribution (2) at- 
tains its maximum value. Note that Vix; Q) > 0. More- 
over, and as a consequence of Eq. (55), we have that 



where 



((17) := ^exp|- l]y(ii;f])|. 



(56) 



(57) 



In this case, Px{x) is a Gibbs distribution with "potential 
energy" function V(x;J7), "temperature" and par- 

tition function C(^)- Clearly, V{x;il.) assigns minimum 
(zero) potential to the states of maximum probability at 
steady-state and infinite potential to the states of zero 
probability. 

We will now assume that, around the thermodynamic 
limit, the potential function V{x; Q) is an analytic func- 
tion of Then, a Taylor series expansion with re- 
spect to around zero (i.e., around the thermody- 
namic limit) results in 

, 1 dV(x;oo) 



Vo{x) + -V^{x) + -- - , 



(58) 



where 



and 



yo(5):=ni;oo) = -^lim lln^>0, 
dV{x; oo) 



V^{x) :- 



As a consequence of Eqs. (56)-(58), and for sufficiently 
large Q, we have that 

Pxix) - exp {- f2Fo(5) - ^i(^)}, 



36 



where the partition function is now given by 
an) = ^exp {- nVoiu) - Vi{u)}. 

u 

Therefore, the stationary probabihty distribution Px(x) 
satisfies a large deviation principle with rate func- 
tion Vo{x) (Touchette, 2009). It turns out that, if 
satisfies the macroscopic equations (44), then 



dVoixjt)) 
dt 



E 



dVM)) dxnjt) 

dXnit) 



dt 



< 0, 



with equality if and only if 



dXnjt) 
dt 



0, for every n S J\f. 



Therefore, the solution x(i) of the macroscopic equa- 
tion (44) produces a downhill motion in the value of the 
potential energy function Vq (which is a Lyapunov func- 
tion for the macroscopic system) until it asymptotically 
reaches a stable stationary state. 

For a given value of O, we can view the multidi- 
mensional surface Vo{x;fl) as a potential energy land- 
scape (Ao et al, 2007; Wang et ai, 2008, 2011; Zhou and 
Qian, 2011). The stable stationary states of Eq. (44) 
correspond to potential wells (basins of attraction) as- 
sociated with the (local or global) minima of Vq, which 
are separated by barriers corresponding to hills (unsta- 
ble states) and saddles (transitional states - states on 
the potential energy surface from which stable states are 
equally accessible). Which minimum will be reached de- 
pends on the initial condition that must be a point in 
the potential well of that minimum. It turns out that, 
once the macroscopic system, given by Eq. (44), reaches 
a stable stationary state, it will stays there forever. As 
a consequence, the macroscopic system will only move 
downhill on the potential energy landscape. 

It is now not difiicult to see that Eqs. (56)-(58) imply 
that 

^ exp I — flVo{x) \ exp < — Vi{x) \ 

lim Px{x) = lim ^ ^-^r ^ 

a^oo a^oo I _ o v^, («) I exp { - Vi («) } 



exp{-Vi(x)} 



if xeg^ 



Euee. cxp{-yi(w)}' 

0, otherwise, 



where ^* is the set of all ground states (global minima) 
of the potential energy landscape Vq (i.e., the set of all 
states at which Vq = 0). As a consequence, only the 
ground states have a non-negligible probability to be ob- 
served as fl becomes large because Px{x) decays expo- 
nentially with fl, for every x ^ Q*. This implies that, 
in the thermodynamic limit, the master equation (5) will 
asymptotically converge almost surely to a ground state 
of the potential energy function Vq, independently of the 



initial state. The particular ground state is chosen with 
probability determined by the values of the potential en- 
ergy function Vi over the ground states of Vq. On the 
other hand, the macroscopic equation (44) will reach a 
minimum of Vq, which may or may not be a ground state 
in 0* depending on the initial condition. If the macro- 
scopic equation has a unique asymptotically stable sta- 
tionary solution that is independent of the initial condi- 
tion, then Vq will have only one (global) minimum. In 
this case, and as we mentioned before, the master equa- 
tion (5) will converge almost surely to the same state in 
the thermodynamic limit. However, if Vq contains more 
than one minimum, then the stationary solution of the 
master equation (5) may be different from the station- 
ary solution predicted by the corresponding macroscopic 
equation (44). As a consequence, 

lim lim px{x;t) lim lim Px{x;t) 

f2— >oo t—>oo t—>oo Q— >oo 

in general. This distinct difference between the station- 
ary behavior of the master equation (left-hand side of 
inequality) and macroscopic equation (right-hand side of 
inequality) is known as Keizer's paradox (Keizer, 1987; 
Qian, 2010; Vellela and Qian, 2007). 

At finite sizes fl, the modes of the stationary proba- 
bility distribution Px{x) (i.e., the most probable states) 
correspond to the minima of the potential energy land- 
scape y(5; Q). If X is a minimum of Vo(x) -|- Q,~^Vi(x), 
then Vq{x) > Vq{x) - VL-^ [Vi{x) - Vi{x)\, for every 
X G W(5')! where >V(5') is a local neighborhood of x 
for which the inequality is satisfied. For large enough O, 
such that 



inax _ max 

x' a;eW(x') 



Vdx)-Vi{x) 
Vo{x') 



< < oo, 



where the first maximum is taken over all minima of 
Vo(^) + n~^Vi{x) which are not ground states of Vq(x), 
we have Vo(x') < Vb(x), for every x € yV(x'), and there- 
fore x' will approximately be a local minimum of the po- 
tential energy landscape Vq. In this case, the peaks of the 
stationary probability distribution Px^) will correspond 
to stationary states of the macroscopic equation (44).^^ 
For this reason, we can refer to the peaks in Px(x) as 
macroscopic modes. 

At smaller values of fJ, the stationary probability dis- 
tribution Pxix) will be given by Eqs. (56) and (57). The 
modes will now depend on the fiuctuation size parame- 
ter and will be determined by the minima of the po- 
tential energy landscape Vix; f2). However, a state that 



Note that the converse of this result is not necessarily true. There 
might be stationary states of the macroscopic equation that, for 
a given value of f2, do not introduce peaks in the stationary 
probability distribution. To sec this, recall that, in the limit 
as r2 — > oo, the only peaks present in the stationary probability 
distribution are the ones associated with the global minima of Vq . 
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minimizes the potential energy function Vq may not nec- 
essarily minimize V, in which case at least some modes of 
the probability distribution Px{x) will not be predicted 
by the corresponding macroscopic equation. Since these 
modes show up at small system sizes, in which apprecia- 
ble stochastic fluctuations may be present in the system 
due to "intrinsic noise," we refer to them as noise-induced 
modes. Recent literature has documented the presence 
of noise-induced modes in biochemical reaction networks 
and their importance in modeling system behavior not 
accounted for by their macroscopic counterparts (Arty- 
omov et ai, 2007, 2009; Bishop and Qian, 2010; Qian, 
2010, 2011; Qian et al., 2009; Zhang et al, 2010c). 

We finally note that, if a Markovian reaction network 
is at a stable state x\ at time to, then it may switch to a 
another stable state at time < t < oo with probabil- 
ity pn{xl]t). However, limo^oo Pn(S2; *) = '^(^2-X(*))i 
where x(t) is the solution of the macroscopic equa- 
tions (42), initialized with x\. Since x\ is a minimum 
of the potential energy function Vq, the macroscopic 
system will be in state x(i) = at time t. Hence, 
limQ_>ooPf!(52j*) = 0. As a consequence, the proba- 
bility of switching from a stable state to another stable 
state tends (in general exponentially) to zero as the sys- 
tem size increases to infinity. At finite system sizes O, 
switching among stable stationary states becomes pos- 
sible, but the probability of switching is very small for 
large O; i.e., switching among stable stationary states 
are rare events (Qian, 2011; Zhou and Qian, 2011). As a 
matter of fact, we can approximate the waiting time for 
switching by an exponential distribution (Aldous, 1982) 
with rate parameter that tends to zero as f2 — >■ oo. 

VIII. MACROSCOPIC (THERMODYNAMIC) BEHAVIOR 

We can view a Markovian reaction network as a ther- 
modynamic system that absorbs energy, produces en- 
tropy, and dissipates heat (Andrieux and Gaspard, 2004, 
2007; Esposito and van den Broeck, 2010a,b; Gaspard, 
2004; Ge, 2009; Ge and Qian, 2010; Oono and Pani- 
coni, 1998; Prigogine, 1978; Qian, 2009, 2010; Ross, 2008; 
Ross and Villaverde, 2010; Santillan and Qian, 2011; 
Schmiedl and Seifert, 2007; Schnakenberg, 1976; Zhang 
et al., 2012). This perspective can provide important in- 
sights into functional properties of Markovian reaction 
networks (such as robustness and stability) and can lead 
to a better understanding of the relationship between 
the mesoscopic (unobservablc) and macroscopic (observ- 
able) behavior of such networks (Gaspard, 2004; Ge et al., 
2012; Han and Wang, 2008; Qian, 2009, 2010; Ross, 2008; 
Ross and Villaverde, 2010; Vellela and Qian, 2009). 

In this section, we consider an irreducible Markovian 
reaction network comprised of M/2 pairs of reversible re- 
actions {2m — 1, 2m), m = 1,2,..., M/2, where 2m — 1 
is the forward reaction and 2m is the corresponding re- 
verse reaction. This does not forbid us to consider ir- 
reversible reactions, since an irreversible reaction can be 
thought of as being reversible with negligible propensity 
in the reverse direction. As we mentioned in Section VI, 



the reaction network is characterized by a unique pop- 
ulation probability distribution px{x;t) that is analytic 
for all t > and converges to a stationary distribution 
Px{x), which does not depend on the initial distribution 
Px{x;0). By following our discussion in Section VII, we 
can define the energy of state x by 



Eix) :-- 



\npx{x), for x G X, 



(59) 



where > is an appropriately chosen size parameter. 

Our discussion in the following is purely mathematical 
in nature and can be applied to any physical or non- 
physical Markovian reaction network. However, direct 
connection to thermodynamics can be made in certain 
physical systems, such as biochemical reaction networks, 
which may exchange matter, work, and heat through a 
well-defined boundary that separates the system with its 
surroundings (Dill and Bromberg, 2011). In this case, we 
must take the size parameter O to be the inverse of fcsT, 
where ks is the Boltzmann constant and T is the system 
temperature. Since the exact value of ^l is not important 
here, we set f2 = 1 for simplicity. 

By viewing a Markovian reaction network as a thermo- 
dynamic system, we can define three fundamental quan- 
tities: the internal energy, entropy, and Helmholtz free 
energy. The internal energy U(t) is the average energy 
of the system at time t over all states, given by 



U{t) ■.= Y,E{x)px{x;t), 



for t > 0,, whereas, the entropy is defined by 



S{t) ■■= - '^Px{x;t)lnpx{x;t), 



(60) 



xex 



for t > 0. Moreover, the Helmholtz free energy is given 

by 



F{t) := f/(i)-5(t) = Vp.(a;;t)ln^ 



xex 



Px{x;t) 



(61) 



for t > 0. The Helmholtz free energy measures the en- 
ergy available in a thermodynamic system to do work 
under constant temperature and volume. Note that 
F{t) coincides with the Kullback-Leibler distance (or 
relative entropy) of the probability distribution px (x; t) 
from the steady-state probability distribution Px (x) [re- 
call Eq. (53)]. Therefore, the Helmholtz free energy pro- 
vides a measure of how far a Markovian reaction net- 
work is from steady-state at time t. Note that F{t) > 
and dF{t)/dt < 0, for every ^ > 0, with equality only 
at steady-state (Cover and Thomas, 1991; Qian, 2009; 
Schnakenberg, 1976). 

A. Balance equations 

Prom Eqs. (5) and (60), we can show the following 
entropy balance equation: 



dS{t) 
dt 



= cj{t) - hit), 



(62) 
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for t > 0, where 

M/2 



(63) 



m=l x^X 



and 



M/2 



m=l xGAf 



-i(a; - S2m-i) 



+Pm(a;;*)in 



7i"2m(a;) 

7''2rTl(a; + S2m-l) 



7r2m-l(a;) 



(64) 



In these equations, 

P+(a;;t) := ■K2m-l{x - S2m-l)Px{x - S2m-i;t) 

p~ (a;; t) := Tr2m{x + S2m-l)Px{x + S2m-i;t) 
- '!r2m-i{x)px{x;t), 

where p'^{x;t) is the net flux of the m-th pair of re- 
versible reactions reaching state x from state x — S2m-i 
and p'^{x;t) is the net flux of the same pair of reac- 
tions reaching state x from state X — 82m (note that 
S2m = -S2m-i)- Moreover, 



-l{x - S2m-l)Pxix - S2m-i;t) 



-^2m{x)Px{x;t) 

'!^2m{x + S2m-l)Px{x + S2m-l]t) 
■K2m-\{x)px{x;t) 



(65) 

are the affinities corresponding to the net fluxes Pm(x; t) 
and p~(a;;f), respectively. Note that 



-A^{x + S2m-i;t), 



and 



dpx (x;t) 
dt 



M/2 



m=l 



for i > 0. Therefore, [p+(x;t) -I- p~(x;t)]dt quantifies 
the change [increase, when p+(x;t) + /9^(x;t) > 0, or 
decrease, when p^{x; t) +p~ (x; t) < 0] in the probability 
mass of the population process within the infinitesimally 
small time interval [t, t + dt) due to the m-th pair of re- 
versible reactions. These changes in probability mass are 
driven by the affinities A^{t) and A^{t), which can be 
viewed as thermodynarnic forces that drive a Markovian 
reaction network away or towards the state of thermody- 
namic equilibrium in which all net fluxes are zero. 

Equation (62) provides an expression for the rate of 
entropy change in a Markovian reaction network. The 
term a{t) quantifies the rate of entropy production. 



whereas, the term h(t) quantifies the rate of entropy 
loss due to heat dissipation. For this reason, (T{t) is 
called the entropy production rate, whereas, h(t) is called 
the heat dissipation rate. Equation (63) shows that 
a{t) is a sum of terms If^J^xex [Pm(^; 0-4m(^; + 
p~ (x; i)^~ (x; t)] , each term quantifying the contribu- 
tion of a pair of reversible reactions to the net rate 
of entropy production. Similarly, Eq. (64) shows that 
h{t) is a simi of terms l/^J^xex {Pm(^; ln[7r2m_i(x - 

S2m-l)/7I"2m(a;)] + p„(x; t) ln[7r2™(x - S2m)/T^2m-l{x)]}, 

each term quantifying the contribution of a pair of re- 
versible reactions to the net rate of heat dissipation. 
Therefore, a reaction with some non-zero net flux must 
produce entropy and dissipate heat. 

By differentiating Eq. (61) with respect to t and by 
using Eqs. (5), (63) and (65), we can show the follow- 
ing balance equations for the Helmholtz free energy and 
internal energy: 



for t > 0, and 



for t > 0, where 



^ M/2 



dFjt) 
dt 



dU{t) 
dt 



fit) - <^(t), 



fit) - hit), 



(66) 



(67) 



fit) -=^1^12 [pt^ix-^t)A+ix) + p;nix;t)A-ix) 

m=l xEX 



with »4+(x) and A^{x) being the affinities of the rn-th 
pair of reversible reactions at steady-state; i.e., A'^{x) 
limt^oo^+(x;t) and A^^ix) := limt_,oo ^;;(x; t). Equa- 
tion (66) quantifles the change in Helmholtz free energy 
due to the Markovian reaction network being away from 
thermodynamic equilibrium at steady-state [quantified 
by the first term on the right-hand-side of Eq. (66)] or 
reduction in Helmholtz free energy due to entropy pro- 
duction [quantified by the second term on the right-hand- 
side of Eq. (66)]. The term f{t) quantifies the rate of 
energy (i.e., power) supplied to the Markovian reaction 
network in order to keep it away from thermodynamic 
equilibrium. For this reason, we refer to /(t) as the 
"motive" power}^ It turns out that < f{t) < a{t), 
for every t > 0.^^ Note that f{t) is a sum of terms 



This quantity is also referred to as "housekeeping" heat rate (Es- 
posito and van den Broeck, 2010a; Ge, 2009; Ge and Qian, 2010; 
Oono and Paniconi, 1998; Schmiedl and Seifert, 2007). However, 
we prefer to call f{t) the "motive" power, since it represents 
the energy flow per unit time required to keep the Markovian 
reaction network away from thermodynamic equilibrium. 
Wc can show the first inequality by using the fact tliat the right- 
hand side of the master equation (5) is zero at steady-state and 
that In a; < x—1, for a; > [see Ge (2009)] . Tiie second inequality 
is due to Eq. (66) and the fact that dF{t)/dt < 0. 



39 



l/2ExeAr [Pm(a;;iM+(x) + p^{x;t)A^{x)], each term 
quantifying the contribution of a pair of reversible reac- 
tions to the net "motive" power. Therefore, a reaction 
with non-zero (forward or reverse) flux and correspond- 
ing non-zero affinity at steady-state will supply motive 
power to the Markovian reaction network. 

Equation (66) shows that reactions in a Markovian re- 
action network can increase the Helmholtz free energy by 
adding "motive" energy to the system, whereas, they can 
reduce the Helmholtz free energy due to entropy produc- 
tion. Moreover, 



a{t) = fit) + 



dF{t) 
dt 



for t > 0, which implies that entropy production comes 

from two soiirccs: from supplying motive power /(t) to 
sustain the reaction network away from thermodynamic 
equilibrium and from a spontaneous change \dF{t)/dt\ 
in Helmholtz free energy due to relaxation towards the 
steady-state (Nicolis and Prigogine, 1977). On the other 
hand, Eq. (67) expresses the first-law of thermodynamics 
(energy conservation): a change AC/ (t) = U{t+dt) — U{t) 
in internal energy within the infinitesimal time inter- 
val [t, t + dt) must equal the amoimt of motive energy 
f{t)dt added to the system minus the dissipated heat 
h{t)dt. Note that a{t) > 0, for every t > 0, with 
equality if and only if A^{t) = A^{t) = 0, for every 
m = 1,2,..., Mjl}^ This is in agreement with the sec- 
ond law of thermodynamics, which postulates that the 
rate of entropy production must always be nonnegative. 
Finally, Eqs. (62) and (66) imply that 



< CT = ft = /, 



(68) 



where a := limj^co o'(^)7 a-nd similarly for h and /. This 
result says that, at steady-state, the amount of motive 
power supplied to the system must be equal to the rate 
of heat dissipation, in agreement with the first law of 
thermodynamics. In addition, the rate of heat dissipation 
must be equal to the rate of entropy production. Finally, 
the steady-state entropy production, heat dissipation and 
motive power must all be nonnegative, in agreement with 
the second law of thermodynamics. 

B. Thermodynamic equilibrium 

A Markovian reaction network roaches thermodynamic 
equilibrium at steady-state if and only if = A^ = 0, 
for every m = 1,2,..., M/2, which is equivalent to the 
following detailed balance equations: 



T^2m{x)Px{x) 



7r2m(a; + S2m-l)Px(^ + «2m-l) = 7r2m-l (a;)Px (^) : 



This is a direct consequence of the fact that (X1 — X2) ln(2;i/rr2) > 
0, for any values of xi and X2, with equaUty if and only if xi = X2- 



for every m = 1, 2, . . . , M/2, x £ X. In this case, f{t) 
I . for every t > 0, which implies that 

'^^'^=-m and ^ = -.(*), 



dt 



dt 



for t > 0. Moreover, Eq. (68) results in a = /i = / = 0, 
which shows that a Markovian reaction network that 
reaches thermodynamic equilibrium at steady-state will 
not produce entropy or dissipate heat. It turns out that 
a Markovian reaction network must be reversible at ther- 
modynamic equilibrium, which means that the stationary 
behavior of the population process will be indistinguish- 
able if the direction of time is reversed. This behavior 
may not be desirable, since many Markovian reaction sys- 
tems (e.g., biochemical reaction networks) are irreversible 
with respect to time. As a matter of fact, entropy pro- 
duction, heat dissipation, and irreversibility with respect 
to time are fundamental properties which are necessary 
for the formation of order in many physical systems (Pri- 
gogine, 1978). Therefore, and in many cases of interest, 
a Markovian reaction network must not reach thermody- 
namic equilibrium in order to be useful. We can make 
sure that this is the case by including nonreversible re- 
actions that transfer mass through the boundary of the 
system with its surroundings, thus breaking detailed bal- 
ance. 

Despite the aforementioned drawbacks, Markovian re- 
action networks that reach thermodynamic equilibrium 
have been extensively used to model population dynam- 
ics. For this type of networks we can use (at least in 
principle) a simple iterative procedure to calculate the 
steady-state probability distribution. This is possible be- 
cause any state x G X can be reached from a given state 
Xo G X through at least one ordered chain of reactions 
(mi, m2, . . . ,mL)- In this case, detailed balance implies 
that (Haken, 1974) 



Px{x) = Px{xo)Y[ 



TTm, [Xo 



1=1 



(69) 



for every x Xq, where mj" is the index of the opposite 
reaction to reaction mi (i.e., = 2m, if m; = 2m — 1, 
and mf = 2m — 1, if m; = 2m). After this procedure 
is completed for all x € A", we can calculate px(^o) in 
Eq. (69) by setting the sum of all probabilities Pxi^) 
equal to 1 —px{xo). 

C. Cycles and affinities 

A useful representation of the state-space A" of a 
Markovian reaction network is by means of a graph G 
whose nodes are the population states x e X and whose 
edges connect pairs of population states (x.x + Sm) 
when 7Tm{x),TTm»{x + Sm) > (Schuakcuberg, 1976). 
Clearly, an edge connecting two states x, x' indicates 
that these states can "reach" each other using a pair 
(m,m*) of reversible reactions. Note that there might 
be a number of distinct edges (corresponding to dif- 
ferent reversible reactions) connecting a given pair of 
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nodes, in which case G is a multi-graph. -'^^ An ordered 
chain (mi, TO2, . . . jUIl) of reactions will produce a path 

{xo,Xo +Smi,- ■■,Xo + Y^fLiSmi) in G of length L, pro- 
vided that each reaction can occur with positive prob- 
ability. In the particular case when "^f-i Smi = 0, the 
reactions (mi, m2, . . . , m^) will produce a cycle in G of 
length L that ends in the same state as the starting state. 
In the following, we will denote by C the set of all cycles 
in G with L > 2." 

It is clear from Eq. (69) that the propensity functions 
of a Markovian reaction network that reaches thermody- 
namic equilibrium must satisfy the following conditions: 



V{C) :=n 



1=1 



1, 



produced by reactions 



over a cycle C = (xo,Xi, . . . ,Xl) 

(mi, m2, ... ,mL), where X; ■= Xq + Y^'i'^i^rn,, ■ These 
are known as Kolmogorov cyclic conditions (Jiang et al., 
2004). Equivalents, 



A{C) ■.= Y^ln- 



1=1 



mi ixi-i)pxixi-i;t) 



= lnP(C) = 0, 



for every t > 0, where A{C) is the net affinity around 
cycle C. In addition to being necessary, the previous 
conditions are also sufficient for a Markovian reaction 
network to reach thermodynamic equilibrium [see Theo- 
rem 2.2.10 in Jiang et al. (2004)]. Therefore, care must 
be taken when dealing with this type of Markovian reac- 
tion networks, since their propensity functions must be 
appropriately constrained. 

Research effort has been recently focused on develop- 
ing techniques for enforcing such conditions in the ther- 
modynamic limit of mass-action systems governed by 
the macroscopic equations (44). In this case, the cor- 
responding constraints are known as Wegscheider condi- 
tions. The proposed methods have been developed for 
performing sensitivity analysis (Zhang et al., 2009) or 
parameter estimation (Colquhoun et at, 2004; Jenkinson 
and Goutsias, 2011; Jenkinson et al., 2010; Liebermeis- 
ter and Klipp, 2006; Yang et al., 2006). However, further 
work is needed to deal with the Kolmogorov- Wegscheider 
conditions in the Markovian setting discussed in this re- 
view. 

The Kolmogorov cyclic conditions are usually highly 
redundant, since a cycle can often be decomposed into 



For example, the reversible reactions Xi X2 and Xi + X3 
X2 + X3 are characterized by the same net stoichiometry and 
will therefore connect the same pair of nodes in G. 
Here we use an equivalence class of cycles over cyclical 

shifts, which implies that a cycle produced by reactions 
(mi,m2, . . . with starting state Xo is equivalent to the cy- 

cle produced by reactions (m2, . . . , mi,, mi) with starting state 
Xo -I- Si. 



smaller cycles. In this case the Kolmogorov cyclic con- 
dition imposed on the larger cycle is implied by the con- 
ditions imposed on the smaller cycles. To address this 
issue, we can derive a minimal set of Kolmogorov cyclic 
conditions that, when satisfied, they imply the remain- 
ing conditions. As a matter of fact, it has been shown 
by Schnakenberg (1976) that the net affinity^(C) of a 
cycle C € C is given by 



K 



A{C) = J2MC)A{Ci), 



(70) 



fe=i 



where ai^{C') is an appropriately defined constant that 
takes integer values and {Cj, C^, . . . , C]^-} is a (non- 
unique) set of cycles, known as fundamental cycles. 

The fundamental cycles and the associated values of a 
can be determined by a simple procedure based on graph 
theory (Schnakenberg, 1976).^^ Given the graph G, we 
can define a maximal tree T{G) of G such that: (i) T(G) 
is a covering subgraph of G [i.e., T(G) shares all ver- 
tices of G and an edge of T(G) must be an edge of G], 
(ii) T{G) is connected, and (iii) T(G) contains no circuit 
(i.e., no cyclic sequence of edges). An edge Ck of G is 
said to be a chord of T(G) whenever it is not an edge of 
r(G). If Tfc(G) is the graph T(G) with the chord Ck in- 
cluded as an edge, then this subgraph of G would include 
exactly one circuit Ck, which is obtained from Tk{G) by 
removing all edges that are not part of the circuit. The 
set {Gi, G2, . . . , Ck} consists of circuits, known as fun- 
damental circuits. Adding an arbitrary orientation to the 
fundamental circuit Gfc defines the fundamental cycle G^ 
used in Eq. (70). 

On the other hand, given a cycle C G C and a funda- 
mental cycle G^, 

a,(G) = ae,(G)ae,(Gt), 

with Cfe being the chord used to produce the circuit Ck 
(and thus G|) and 

i, if cycle D contains i copies of an 
edge e in the same orientation 
as in graph G 

Ge{D) := ^ —i, if cycle D contains i copies of an 
edge e in the opposite orientation 
as in graph G 

0, if cycle D does not contain edge e, 

where the graph G is assigned an orientation to its edges 
in the direction of the forward reaction (i.e., the reaction 
that corresponds to an odd value of m). 



18 ■\Yc refer the reader to the excellent book by Diestel (1997) for 
an introduction to graph theory. 
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Equation (70) shows that the affinity aronnd a cycle C 
can be written as a Unear combination of the affinities 
around the fundamental cycles C'l- This result demon- 
strates that a necessary and sufficient condition for a sys- 
tem to reach thermodynamic equilibrium is that the affin- 
ity of each fundamental cycle must be zero. This provides 
a reduced and more manageable set of conditions than 
the Kolmogorov cyclic conditions over all possible cycles. 

The not affinity of a Markovian reaction network that 
does not reach thermodynamic equilibrium must be non- 
zero over at least one fundamental cycle. This affin- 
ity quantifies the net thermodynamic force applied to 
the network due to its interaction with the surroundings 
(e.g., due to mass flow through the system boundary); 
see Andrieux and Gaspard (2004, 2007) and Schnaken- 
berg (1976). In many Markovian reaction networks, 
such as those with propensities that follow the mass- 
action rate law, there is a number of global affinities Aq, 
q = 1,2,. . . ,Q, which describe the macroscopic coupling 
of the system to its surroundings, such that the affin- 
ity of each fundamental cycle equals Aq, for some q. 
If the global affinities are known, then the relationships 
A{Cl) = Aq, along with Eq. (70), constitute a more 
general version of the Kolmogorov cyclic conditions that 
must be enforced on the propensity functions so that the 
Markovian reaction network does not reach thermody- 
namic equilibrium. 

D. Example: Neural dynamics 

We now consider a special case of the neural network 
model discussed in Section III-F, which allows us to nu- 
merically compute the dynamics of the joint population 
probability distribution and proceed with illustrating a 
number of thermodynamic properties of this model. We 
will assume that the L neurons in the network can be 
divided into an equal number of L/2 excitatory and L/2 
inhibitory neurons. Moreover, we will assume that all 
neurons synapse to all other neurons with excitatory and 
inhibitory weights i^e ^ ^ and I'j < 0, respectively. Fi- 
nally, wo will consider the case in which the neurons are 
characterized by the same decay rate 7, whereas, their 
external inputs take the same value h. Under these as- 
sumptions, it may not be of interest to track the state 
of individual neurons, since the excitatory or inhibitory 
neurons are identical to each other. Instead, it will be 
more appropriate to track the dynamics of the net num- 
ber 

Ait) :=Y,{t)+Y2{t), 
of active excitatory and inhibitory neurons, where 
Y^it) := J2 ^2i{t) and Y^it) := ^ X^iit), 

with £ and I being the set of excitatory and inhibitory 
neurons, respectively. 

For convenience, and without loss of generality, we 
can take £ = {1,2,..., L/2} and I = {L/2 + 1, L/2 + 



2, . . . , L}. It turns out that Yi{t) and Y2{t) can be mod- 
eled by a simple Markovian reaction network comprised 
of two species Yi and Y2 that denote active excitatory and 
inhibitory neurons, respectively, which interact through 
the following reactions: 

Yi+Y2^ 2Yi + Y2 

Fi + 1^2 ^ Fi + 2^2 
Y2 0. 

These reactions correspond to the activa- 
tion/deactivation of an excitatory neuron (first and 
second reactions) and the activation/dcactivation of an 
excitatory neuron (third and fourth reactions). Their 
propensity functions are given by [recall Eqs. (12) 
and (13)] 

7ri(3/) = (L/2 - yiMy) > 0] tanh[,^(y)] 
7r2(3/) = 72/1 

7r3(y) = {L/2 - y2My) > 0] tanh[.^(j,)] 

7I"4(3/) = 72/2, 

respectively, where y = [y\ 2/2]^, with yi and 1/2 taking 

values in {0, 1, ... , L/2}, whereas, [a > 0] is the Iverson 
bracket and (j) is the synaptic input to each neuron, given 
by [recall Eq. (11)] 

(l>(y) = vevi + vm + h. 

The propensity functions of the associated reverse reac- 
tions are taken to be zero. 

Despite its simplified nature, the previous model has 
been shown by Benayoim et al. (2010) to be very effec- 
tive for predicting experimentally observed, in vitro and 
in vivo, neural behavior, known as avalanches. This be- 
havior is characterized by irregular and isolated bursts 
of neural activity during which many neurons fire simul- 
taneously. In the following, we use the thermodynamic 
principles discussed in this section to explore this inter- 
esting behavior. For ease of computational analysis, we 
consider a moderately sized neural network comprised of 
L = 100 neurons. This allows us to numerically compute 
the solution of the imdcrlying master equation using the 
KSA method discussed in Section IV-B. We adopt pa- 
rameter values used in Benayoun et al. (2010) and set 
7 = O.lms"^ and h = 0.001, whereas, we chose values 
for the synaptic weights i^e and vj so that their sum 
i^E + is kept fixed to a value of 0.004. Finally, we as- 
sume that all neurons are initially at rest, in which case, 
yi(0) = F2(0)=0. 

It has been shown by Benayoun et al. (2010) that 
when the sum ve + vi is kept fixed, then the difference 
5v := Ve — vi controls avalanching (see also Fig. 10). 
More particularly, as Ev increases, the network transi- 
tions from asynchronous to synchronous neural firings 
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FIG. 8 Dynamic evolutions of the net number of active 
neurons in the neural network example considered in Sec- 
tion VIII-D drawn from the underlying master equation using 
exact sampling. The blue trajectory indicates that neurons 
fire asynchronously, whereas, the irregular and isolated bursts 
of net neural activity observed in the red trajectory indicate 
that neurons fire synchronously resulting in avalanching. 

that lead to avalanching. We iUustrate this behavior in 
Fig. 8, which depicts two trajectories of the net number 
A{t) of active neurons obtained by samphng the master 
equation using exact sampling. The blue trajectory has 
been obtained by setting ve = 0.034 and i/j — —0.00062, 
whereas, the red trajectory has been obtained by set- 
ting ve — 0.140 and i^i = —0.136. Clearly, the blue 
trajectory indicates that neurons fire asynchronously (in 
this case, 6v = 0.00276), whereas, the red trajectory ex- 
hibits avalanching with neurons firing synchronously, re- 
sulting in irregular and isolated bursts of activity and 
thus avalanching (in this case, Siy = 0.276, which is 100 
times larger than the previous value). 

Most biological systems of interest reach a state of 
homeostasis, wherein the system is maintained at a given 
stable operating point. Mathematically, we can describe 
this point by a stable steady state. From a thermody- 
namic perspective however such a system must operate 
away from thermodynamic equilibrium, since living or- 
ganisms require transfer of energy and mass with their 
surroundings in order to consume nutrients and excrete 
waste. To achieve this, nonzero motive power must be 
supplied to the system at steady-state, which implies that 
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FIG. 9 Dynamic evolutions of internal energy, entropy, 
Helmholtz free energy, entropy production rate, heat dissi- 
pation rate, and supplied motive power in the neural network 
example considered in Section VIII-D for the case of asyn- 
chronous (blue lines) and synchronous (red lines) neural fir- 
ings leading to avalanching. 

an equal amount must be dissipated to the surroundings 
in the form of heat, by virtue of the first law of thermo- 
dynamics. As a consequence, the state of homeostasis at 
which biological systems operate is often referred to as 
the non- equilibrium steady-state (NESS). 

Naturally, the neural network discussed in this exam- 
ple must also operate at a NESS [e.g., see Stewart and 
Plenz (2008)]. Fig. 9 shows clearly that this is indeed 
the case. This figure depicts the dynamic evolutions of 
internal energy, entropy, Helmholtz free energy, entropy 
production rate, heat dissipation rate, and supplied mo- 
tive power, for the two cases of asynchronous (blue lines) 
and synchronous neural firings (red lines). Note that 
all thermodynamic quantities reach stationarity, with the 
entropy production rate, heat dissipation rate, and sup- 
plied motive power converging to the same value in each 
case, in agreement with the first and second law of ther- 
modynamics. The fact that this value is nonzero in both 
cases shows that the system operates away from ther- 
modynamic equilibrium at steady-state regardless of the 
type of neural firings involved. 

The results depicted in Fig. 9 show that the system 
entropy is in general smaller when neurons fire syn- 
chronously than when they fire asynchronously. This 
indicates an expected degree of predictability in neural 
activity when neurons fire synchronously. On the other 
hand, the initial Helmholtz free energy associated with 
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asynchronous neural firings is almost an order of mag- 
nitude larger than that associated with synchronous fir- 
ings. In both cases however the Helmholtz free energy 
becomes zero at steady-state, as expected. This indicates 
that, under constant temperature and volume, apprecia- 
ble more work must be done when the neurons fire asyn- 
chronously to reach steady-state than when the neurons 
fire synchronously. It is therefore expected that when 
neurons fire synchronously, the system will reach steady- 
state much faster than when neurons fire asynchronously. 
This is clearly verified by the results depicted in Fig. 9. 
We may therefore postulate that synchronous neural fir- 
ing is, among other things, necessary for a neural network 
to quickly reach a state of homeostasis [see also Stewart 
and Plenz (2008)]. 

Fig. 9 also reveals that the stationary value of the sup- 
plied motive power (as well as the stationary values of 
the entropy production and heat dissipation rates) is ap- 
preciably larger in the case of synchronous neural firings 
than asynchronous firings. This difference is well pre- 
dicted by the theory of dissipative structures (Nicolis and 
Prigogine, 1977) according to which self-organization of 
a system to an ordered internal state requires that the 
system is sufficiently driven by external sources and dis- 
sipates appreciable heat to its surroundings. It is clear 
from Fig. 9 that appreciable amounts of supplied mo- 
tive power and heat dissipation is required to achieve 
avalanching behavior, which leads to an ordered station- 
ary state, quantified by a lower entropy. We may there- 
fore conclude that the emergence of avalanches in a neu- 
ral network is a consequence of externally driven self- 
organization accompanied by appreciable heat dissipa- 
tion. This conclusion is further confirmed by Fig. 10, 
which depicts the (average) rate of avalanche formation 
(number of avalanches per unit time, ) the supplied mo- 
tive power (or heat dissipation), and the system entropy 
at steady-state. To calculate the number of avalanches 
present in a given trajectory of net neural activity, we 
assume that an avalanche occurs within a time window 
[t, t + r) whenever the following three conditions are sat- 
isfied: (z) A{t') > 0, for aU t' e [t,t + t); i.e., there 
is neuronal activity during the time interval [t,t + r), 
{ii) there exist some e > such that A{t') = 0, for all 
t' G [t — e, t); i.e., there is no neural activity immediately 
before time t, and (iii) A{t + t) = 0; i.e., there is no 
neural activity at time t + r. 

We conclude by noting that avalanching can be un- 
derstood at steady-state through the energy landscape 
E{y) [see Eq. (59)]. As a matter of fact, the stationary 
dynamics of neural activity may be thought of as a ran- 
dom walk on this landscape, where the most likely steps 
follow a path from higher to lower energy states (i.e., 
the preference is to move downhill), with occasional (low 
probability) jumps from lower to higher energy states. In 
Fig. 11(a), we depict the energy landscape when neurons 
fire asynchronously, whereas, in Fig. 11(b) we depict the 
energy landscape when neurons fire synchronously, thus 
leading to avalanching. One can see that the ground 
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FIG. 10 Average rate of avalanche formation (red curve) in 
the neural network example considered in Section VIII-D (cal- 
culated from 1,000 trajectories of net neural activity during a 
period of 1,000 msec, drawn from the master equation using 
exact Monte Carlo sampling), supplied motive power (or heat 
dissipation) at steady-state (green curve), and system entropy 
at steady-state (blue curve) as a function of the difference 5v 
between the excitatory and inhibitory weights. As expected, 
increasing the supply of motive power results in increasing 
the rate of avalanche formation and decreasing the system 
entropy. 

state (global minimum) of the energy landscape depicted 
in Fig. 11(a) occurs away from the origin at which neu- 
ral activity is zero. As a consequence, the system spends 
most time in Gaussian-like fiuctuations about the ground 
state. The system can jump out of the energy well sur- 
rounding the ground state and reach the origin, but with 
very small probability, due to its large width. Therefore, 
avalanche formation in this system is a rare event. This 
behavior is in agreement with a recent finding that neu- 
ral networks may simultaneously support synchronous 
and asynchronous dynamics, switching between these two 
modes of operation spontaneously (De Ville and Peskin, 
2008). 

On the other hand, the energy landscape of the system 
with synchronous neural firings depicted in Fig. 11(b) 
contains a valley along the line Yi — Y2 of equal excita- 
tory and inhibitory activities, which slopes down to the 
ground state that is now located at the origin. Thus, from 
any point on this landscape, the most likely trajectory 
roles downhill until it reaches the origin. From this point, 
the dynamics of excitatory and inhibitory neural activ- 
ities may again reach the valley by randomly jumping 
away in an uphill motion that overcomes the steep and 
narrow energy well surrounding the origin. This mecha- 
nism results in avalanching, which can be thought of as 
a random sequence of zero (or almost zero) net neural 
activity at points in the well surrounding the origin fol- 
lowed by nonzero net neural activity at points proximal 
to the valley. 
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FIG. 11 Energy landscape of the neural network example considered in Section VIII-D when: (a) neurons fire asynchronously, 
and (b) neurons fire synchronously resulting in avalanching. 



IX. OUTLOOK 

The study of Markov processes on complex reaction 
networks is an area of research that has been evolving 
for more than half a century. Its applicability to many 
scientific and engineering disciplines has led to parallel 
and often independent developments, which have recently 
reached a critical mass due to unprecedented advance- 
ments in modern experimental procedures and computa- 
tional capabilities. This can be best illustrated by the 
sharp increase in the number of articles published on the 
subject during the last ten years. For this reason, we 
believe that the present review is both appropriate and 
timely. Our objective has been to provide a coherent 
exposure to what has been done so far and illustrate 
key methods with simple examples. In the following, 
we conclude this review with a brief discussion of some 
outstanding issues and problems important to the field. 
Some old problems require better solutions, whereas, de- 
veloping novel methods for the analysis of Markovian re- 
action networks can lead to new and important results. 
In some cases, the work that needs to be done is at least 
as challenging and rewarding as the work done so far. 

A. Solving the master equation 

The tremendous flexibility and generality of Markovian 
reaction networks make them an excellent mathematical 
framework for studying stochastic processes on complex 
networks. The coherency of a single framework means 
that tools and discoveries made in one field may be read- 
ily ported to distant applications. In fact, it has been 
argued by Cook et al. (2009) that Markovian reaction 
networks with mass action propensities can perform Tur- 
ing universal computations with arbitrarily small error, 
which becomes zero at the thermodynamic limit (Mag- 
nasco, 1997). This strength however turns out to be 



one of the most profound weaknesses of Markovian re- 
action networks: there will be no single analytical or 
even computational method capable of calculating the 
exact solution of the underlying master equation in com- 
plete generality using finite resources. As a consequence, 
the development of accurate and computationally feasi- 
ble techniques for studying the dynamic behavior of large 
nonlinear Markovian reaction networks is still the most 
important and challenging problem in this field of re- 
search. 

One way to deal with this problem is to focus on spe- 
cialized structures that may be present in certain reac- 
tion networks and, by exploiting these structures, develop 
rigorous approximation techniques tailored to the specific 
application at hand. A good example of such an approach 
is the IE method discussed in Section IV-B. This method 
calculates the exact solution of the master equation, up 
to a desired precision, by exploiting the structure of the 
master equation that governs the DA process. More- 
over, it uses the fact that, in some reaction networks, the 
sample space associated with the DA process is bounded 
whereas its cardinality is not appreciably larger than the 
cardinality of the sample space associated with the pop- 
ulation process (Jenkinson and Goutsias, 2012). 

Another possibility is to move away from estimating 
the joint probability distributions of the DA and pop- 
ulation processes and focus on estimating computation- 
ally more tractable marginal probability distributions or 
statistical summaries. Moment closure schemes, in con- 
junction with MaxEnt methods, seem to be particularly 
suited in this case. However, much work is needed for de- 
veloping appropriate closure schemes and evaluating the 
errors introduced by the resulting approximations, as well 
as for designing computationally efficient MaxEnt meth- 
ods. On the other hand, embedding Monte Carlo steps 
within a computationally efficient method for solving the 
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master equation can prove very useful for guiding cor- 
rect implementation. For example, employing a second- 
order moment closure scheme will not be appropriate if 
a small number of Monte Carlo samples drawn from the 
master equation reveals a bistable system behavior. The 
use of hybrid techniques, capable of drawing on relative 
strengths to mitigate weaknesses of the tools discussed in 
Sections IV & V, will likely pave the way to more robust 
solution methodologies, especially for stiff reaction net- 
works. Finally, to deal with the large networks present 
in many applications, it is necessary to focus on compu- 
tational efficiency and, in particular, on algorithms that 
can exploit the highly parallel structure of modern high 
performance computing platforms. 

B. Thermodynamic analysis 

Statistical thermodynamics can be effectively used to 
compactly describe the macroscoi)ic b(4iavior of a given 
stochastic system by a small number of variables. For 
this reason, it has been successfully applied in many 
fields of science and engineering. For example, a con- 
tainer full of gas molecules can be exhaustively de- 
scribed by a set of high-dimensional Hamiltonian equa- 
tions governing the position and momentiim of every sin- 
gle molecule. However, a small number of statistical ther- 
modynamic summaries, such as pressure, temperature, 
entropy, Helmholtz free energy etc., can be used to pro- 
vide a more lucid and computationally tractable descrip- 
tion of the system. Compact descriptions of system dy- 
namics are also possible in the case of Markovian reaction 
networks whose analysis, based on statistical thermody- 
namics, may be the only amenable method of dealing 
with large nonlinear networks. However, development of 
such analysis methods are still at their infancy and wide 
open for future exploration. 

A promising line of inquiry seems to be the potential 
energy landscape perspective discussed in Section VH, 
which serves as the starting point for developing the ther- 
modynamic analysis tools discussed in Section VHI. Fu- 
ture efforts must focus on designing accurate and efficient 
methods for estimating the potential energy landscape 
of Markovian reaction networks and detecting noise- 
induced modes. The ability to predict noise- induced 
modes not present as fixed points in the macroscopic 
equations is an important and challenging task. Since 
these modes of operation are prominent when appreciable 
stochastic fluctuations are present and, often, arc perched 
near bifurcation points predicted by the corresponding 
macroscopic equations (Samoilov et ai, 2005; Turcotte 
et al., 2008), it is likely that their study will lead in- 
vestigators to focus on the interface between stochastic 
processes and bifurcation theory for nonlinear ODEs. 

Aside from the possibility of representing a reaction 
network using a small number of variables, conditions on 
the underlying propensity functions imposed by the laws 
of thermodynamics have shown to be extremely valuable 
for the modeling and analysis of such networks. For ex- 
ample, the use of thermodynamic constraints can allevi- 



ate some burden associated with parameter estimation 
imposed by the curse of dimensionality and can lead to 
lower computational complexity, better estimation per- 
formance, and reduced data overfitting (Colquhoun et al, 
2004; Jenkinson and Goutsias, 2011; Jenkinson et al, 
2010; Liebermeister and Klipp, 2006; Yang et al., 2006). 
Moreover, it can result in physically realizable network 
models consistent with the fundamental laws of thermo- 
dynamics. The work done so far on this important sub- 
ject has focused on reaction networks with deterministic 
dynamics. Therefore, extending this line of research to 
Markovian reaction networks is an exciting prospect with 
potentially fundamental consequences. 

C. Sensitivity analysis 

Often, the main focus of analysis of the dynamic be- 
havior of a reaction network is a response function that 
encapsulates some important system characteristics. In 
epidemiology, for example, one may not care so much 
about the specific details of the population dynamics, but 
would rather focus on the total number of individiials in- 
fected by a disease over a given period of time. Another 
example would be the case of cell signaling, where the 
detailed interactions of a signaling pathway are not as 
important as the total amount of a protein produced at 
the "output" of the pathway. Sensitivity analysis is a 
quantitative approach designed to investigate how vari- 
ations in the parameters of a reaction network (e.g., in 
the specific probability rate constants associated with the 
propensity functions of a mass action system) affect a re- 
sponse function of interest (Heinrich and Schuster, 1996; 
Saltelh et al, 2008, 2005; Varma et at, 1999). 

Many physical and man-made reaction networks are 
designed to be robust to random fluctuations (or even 
failures) in system components. Although robustness is 
a highly desirable property, it results in a small number 
of parameters having a disproportionately large influence 
on the system response. As a consequence, a robust reac- 
tion network can be quite vulnerable to targeted attacks 
on influential components, which can be a blessing or a 
curse, depending on the particular situation at hand. For 
example, development of new drugs may greatly beneflt 
from this property since, to reduce or even eliminate the 
effects of a disease caused by deregulation of key system 
responses, it may be sufficient to design a drug that only 
inhibits influential reactions that shape these responses. 
On the other hand, targeted attacks on national infras- 
tructure by hackers or terrorists may produce large scale 
disruptions with devastating results. 

The objective of sensitivity analysis is to determine 
those factors in a reaction network that produce no no- 
ticeable variations in system response and identify those 
factors that are most influential in shaping that response. 
Although this is a powerful analysis technique with im- 
portant practical consequences, it comes with a large 
computational cost, even in the case of reaction net- 
works with deterministic dynamics (Zhang et al, 2009; 
Zhang and Goutsias, 2010, 2011). For this reason, the 
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development of practical methods for sensitivity analy- 
sis of Markovian reaction networks is still at their in- 
fancy (Dandach and Khammash, 2010; Gunawan et ai, 
2005; Kim et al, 2007; Kim and Sauro, 2010; Komorowski 
et al., 2011; Plyasunov and Arkin, 2007; Rathinam et al., 
2010; Warren and Allen, 2012). 

In the stochastic context, sensitivity analysis involves 
computing the solution of the master equation using dif- 
fercmt parameter values. As a consequence, the devel- 
opment of efficient solution methods that can be imple- 
mented on parallel computer architectures, paired with 
novel sensitivity estimators, will ensure the feasibility of 
this type of analysis. We should also note that it has been 
recently demonstrated by Zhang et al. (2009) and Zhang 
and Goutsias (2011) that, at least for the case of physical 
reaction networks with deterministic dynamics, sensitiv- 
ity analysis methods must be in agreement with under- 
lying thermodynamic constraints. As a consequence, de- 
veloping accurate, computationally efficient, and thermo- 
dynamically consistent sensitivity analysis methods for 
Markovian reaction networks is an important research 
activity with significant benefits. 

D. Statistical inference 

In general, there are two fundamentally different types 
of parameters associated with a Markovian reaction net- 
work model: the stoichiometric coefficients Vnm and f'nm 
that determine the structure of the network, and the ki- 
netic parameters that determine the non-structural por- 
tion of the propensity functions. Some parameter val- 
ues can be deduced experimentally or by means of ap- 
propriate theoretical and sometimes heuristic arguments. 
Most parameters however must be estimated from avail- 
able data using statistical inference techniques. Since 
the predictive power of a given model is fundamentally 
constrained by the accuracy of its parameterization, in- 
ferring the unknown parameter values in a Markovian 
reaction network is a problem of paramount interest and 
practical importance. Although this problem has been 
extensively studied for reaction networks with determin- 
istic dynamics (Crampin et al., 2004; Maria, 2004; Moles 
et al., 2003), the statistical inference of Markovian reac- 
tion networks is largely an open research problem. This 
problem has been recently investigated by Golightly and 
Wilkinson (2005, 2006), Reinker et al. (2006), Boys et al. 
(2008), Komorowski et al. (2009), Poovathingal and Gu- 
nawan (2010), Wang et al. (2010d), and Daigle Jr. et al. 
(2012), but the resulting algorithms do not adequately 
address important issues, such as curse of dimensional- 
ity, thermodynamic consistency, and computational effi- 
ciency. These methods have been primarily designed for 
biochemical reaction networks, but can be easily adopted 
in other applications with little or no effort. 

In most approaches to statistical inference, it is quite 
common to assume known structural parameters and pro- 
ceed with estimating the kinetic parameters using noisy 
and sparse measurements of system dynamics. This 
problem, known as model calibration, is much easier than 



the problem of estimating the structural parameters, 
which is often referred to as model selection. 

The two most difficult issues associated with model 
calibration is the curse of dimensionality and the use 
of non-convex cost functions which complicate numer- 
ical optimization. Curse of dimensionality refers to a 
very fast growth of the volume of the parameter space in 
terms of the number of unknown parameters to be esti- 
mated. As a consequence, the problem of finding the 
"best" parameter values becomes increasingly difficult 
as the number of unknown parameters increases. This 
is further exacerbated by the non-convex optimization 
problem of finding these values, which is computation- 
ally very difficult to solve in most cases of interest (Spall, 
2003). Therefore, the development of statistical tech- 
niques for accurate and computationally efficient model 
calibration of Markovian reaction networks is a very chal- 
lenging problem. Possible ways to attack this problem is 
to effectively reduce the number of parameters that must 
be estimated by incorporating appropriate constraints 
[e.g., constraints imposed by the fundamental laws of 
thermodynamics (Colquhoun et al., 2004; Jenkinson and 
Goutsias, 2011; Jenkinson et al., 2010; Liebermeister and 
Klipp, 2006; Yang et al, 2006)] and to identify a smaller 
set of "influential" parameters whose values must be es- 
timated with sufficient precision [e.g., by employing a 
sensitivity analysis approach (Jenkinson et al., 2010)]. 
This reduction in dimensionality must be combined with 
fast algorithms for solving the master equation, with effi- 
cient optimization methods, and appropriately designed 
experimental protocols for collecting data with high in- 
formation content about the values of the unknown pa- 
rameters (Jenkinson and Goutsias, 2010). 

In general, model selection is a more difficult prob- 
lem. Solving this problem will require the development 
of novel hypothesis testing approaches for comparing be- 
tween two competing network models (e.g., an originally 
proposed signaling network and another network ob- 
tained by adding new reactions) in a rigorous statistical 
fashion. This approach however requires that both mod- 
els are calibrated before compared to each other (e.g., by 
a likelihood ratio test), which substantially adds to the 
difficulty of the problem. Another major issue is that 
more complex models are expected to be more capable 
of closely matching experimental data, but these models 
may result in undesirable overfitting. It is therefore nec- 
essary to develop methods that appropriately penalize 
model complexity so that the chosen "optimal" model is 
the most parsimonious model capable of adequately ex- 
plaining available data. Finally, all of this must be done 
while taking into account possible constraints imposed 
on the structural and kinetic parameters of the network 
(e.g., by prior knowledge on feasible structural parameter 
values and by the fundamental laws of thermodynamics). 

E. Adaptive IVIarlcovian reaction networlcs 

An important aspect of many real-life interaction net- 
works is that the network topology, quantified by the sto- 
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ichiomctric cocfRcients and v'^m^ intricately cou- 
pled with the network dynamics. As we mentioned in 
Section II-B-3, for all reaction networks enconntered in 
practice, the propensity functions depend on the stoi- 
chiometric coefficients Vnm and, therefore, the topolog- 
ical structure of a reaction network directly affects its 
dynamics. However, the opposite is also true: the dy- 
namics can influence the underlying network topology. 
In opinion formation, for example, individuals form their 
beliefs based on interactions determined by the under- 
lined social network, whereas, individuals with similar 
beliefs tend to eventually influence each other. In these 
cases, modeling dynamics on complex networks by as- 
suming fixed topology is legitimate only when the time 
scale of interest is sufficiently smaller than the time scale 
of change in network topology. 

Recently, several articles have appeared in the lit- 
erature introducing adaptive networks that take into 
account the interplay between network topology and 
dynamics (Ehrhardt et al, 2006; Gross and Blasius, 
2008; Gross et al., 2006; Holme and Newman, 2006; 
Ren et al, 2010; Yuan and Zhou, 2011; Zhang et al, 
2010b; Zhou and Kurths. 2006). These preliminary works 
clearly demonstrate that a number of intriguing prop- 
erties emerge, not previously observed in nonadaptive 
networks: formation of complex topologies, spontaneous 
emergence of modular organization, more complex dy- 
namics than the ones observed in nonadaptive models, 
and self-organization towards a highly robust critical be- 
havior characterized by power-law distributions. 

Understanding the coupling between Markovian dy- 
namics and network topology is a fascinating subject of 
research that can eventually lead to surprising new the- 
oretical results and novel computational methods for the 
modeling and analysis of networks with significant im- 
pact in many fields. For example, adaptive Markovian 
reaction networks may potentially lead to a better un- 
derstanding of how cell biology at the molecular level has 
evolved towards certain types of reaction network archi- 
tectures and motifs, characterized by a high degree of 
modularity and robustness. These type of networks can 
be represented by master equations with timc-dcpcndcnt 
stoichiometric coefficients, whose values are updated by 
the system state following well-defined rules. We foresee 
a rapid growth in developing such models that can lead 
to many fascinating and novel results to come. 

X. CONCLUSION 

The recent burst in the amount of research efforts for 
studying random processes on complex interaction net- 
works has been driven primarily by impressive advances 
in experimental techniques for measuring these processes 
and by a clear understanding that stochasticity plays a 
fundamental role in shaping the transient and steady- 
state dynamic behavior of real-life networks. This ef- 
fort has reinforced the fact that the theory of Markovian 
reaction networks is at the foundation of most model- 
ing and analysis techniques for studying stochastic dy- 



namical systems on networks with applications in diverse 
fields, such as chemistry, biology, sociology, epidemiology, 
pharmacology, theoretical neuroscience, and engineering. 
The main goal of this review was to summarize, in a 
systematic fashion, what is known so far about this ex- 
citing field and briefly discuss some open problems that 
need to be solved and new methodologies that must be 
developed. We believe that a major effort should be fo- 
cused on introducing new concepts and ideas that might 
lead to novel and practical tools for the modeling and 
analysis of stochastic dynamics on interaction networks 
encoimtercd in practice. We hope that this article will 
be used as a reference by network scientists across dif- 
ferent scientific disciplines and help to catalyze exciting 
new developments in these fields. 
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